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1  Summary 

This  document  summarizes  the  results  of  a  United  States  Air  Force  Small  Business  Innovation 
Research  (SBIR)  project  entitled,  “Multiparticle  Impact  Damping  Design  Methodology  for 
Extreme  Environments.”  The  intent  of  this  project  was  to  develop  and  validate  a  method  for 
designing  particle  damping  into  complex  structures. 

Multiparticle  impact  damping  (MPID),  or  particle  damping,  is  a  derivative  of  impact  damping 
where  multiple  auxiliary  masses  of  small  size  are  placed  in  a  cavity  attached  to  the  vibrating 
structure.  Particle  damping  can  perform  at  elevated  temperatures  where  most  other  forms  of 
passive  damping  cannot.  Studies  conducted  over  recent  years  have  demonstrated  the 
effectiveness  and  potential  application  of  particle  dampers,  and  have  shown  that  particle 
dampers  are  highly  nonlinear  dampers  whose  energy  dissipation,  or  damping,  is  derived  from  a 
combination  of  loss  mechanisms.  The  relative  effectiveness  of  these  mechanisms  changes 
based  on  various  system  parameters. 

A  multi-tiered  approach  was  used  to  develop  the  design  methodology: 

•  Perform  experimental  studies  to  characterize  particle  damping  on  vibrating  structures. 

•  Examine  changes  in  particle  damping  effectiveness  due  to  environment. 

•  Develop  an  analytical  modeling  capability  to  determine  how  parameters  in  a  particle 
damper  contribute  damping  to  a  vibrating  structure. 

Due  to  the  complex  interactions  of  the  loss  mechanisms  in  a  particle  damper  and  the  large 
number  of  parameters  affecting  the  damper  performance,  it  is  extremely  difficult  to  explicitly 
define  a  particle  damper  configuration  for  a  particular  application.  However,  based  on  the 
damper  behavior  observed  in  experimental  testing  and  analytical  simulations,  design  guidelines 
have  been  established. 
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2  Introduction 

Vibration  mitigation  is  often  required  in  aerospace  structures,  including  engine  subsystems.  Yet 
established  damping  methods  are  not  effective  in  high-temperature  environments.  MPID  is  a 
promising  technology  that  can  be  effective  over  a  wide  temperature  range.  Energy  is  dissipated 
as  particles  within  a  damping  system  impact  one  another  and  the  walls  of  their  enclosure.  The 
actual  design  of  this  type  of  damping  is  complex  and  this  research  focused  on  development  of  a 
design  methodology.  Analytical  modeling  capability  is  central  to  the  methodology,  but  a  purely 
analytical  approach  is  not  feasible. 

MPDD  consists  of  many  small  metallic  or  ceramic  particles  contained  in  a  cavity  that,  when 
excited  by  vibrations,  cause  the  base  structural  motion  to  be  damped.  The  impact  of  particles  on 
each  other  and  on  the  cavity  walls,  friction  between  particles,  and  friction  between  particles  and 
the  cavity  walls  causes  energy  dissipation  which  reduces  the  amplitude  of  the  base  structure 
vibration.  Since  the  particles  can  be  metallic  or  ceramic,  they  can  be  used  at  high  temperatures, 
and  since  they  can  be  designed  to  be  temperature  insensitive,  the  damping  mechanisms  are 
fairly  temperature  insensitive  over  wide  temperature  ranges. 

Research  has  indicated  that  MPID  technology  may  offer  a  robust  solution  to  vibration  problems 
under  conditions  where  other  damping  solutions  are  not  viable.  Viscoelastic  materials  and 
vitreous  enamels,  while  presenting  well-developed  and  effective  solutions  to  problematic 
vibration,  do  not  offer  practical  damping  solutions  at  either  high  temperatures  or  over  extreme 
temperature  ranges.  Friction  dampers  may  have  limited  life  expectancies  due  to  the  surface 
effects  changing.  Viscous,  magnetic,  and  passive  piezoelectric  damping  cannot  perform  at  high 
temperatures.  The  loss  mechanism  for  powder  damping  may  be  a  damping  mechanism  that  will 
provide  damping  over  high  temperatures  or  over  extreme  temperature  ranges.  However,  this 
mechanism  is  a  viscous  flow  mechanism,  with  the  powder  under  high  compressive  loads.  This 
means  that  one  needs  either  a  device  that  can  replenish  the  powder  as  it  flows  or  the  vibration 
amplitude  must  be  sufficient  to  overcome  the  compressive  loads  to  flow  back.  This  leaves 
multiparticle  damping  as  the  best  loss  mechanism  for  extreme  temperature  environments. 

Previous  work  in  impact  damping  forms  the  basis  for  new  developments  in  MPID.  Impact  of 
one  element  with  another  can  be  modeled  reasonably  effectively,  and  various  impact  dampers 
have  been  designed  and  built  over  time.  Of  the  several  efforts  to  develop  and  deploy  multiple 
particle  dampers,  some  have  been  successful,  but  there  is  no  clear  underlying  design 
methodology  that  enables  rapid  accurate  design  for  new  applications. 

Modeling  of  particle  damping  is  difficult  because  of  the  inherent  nonlinearities  in  the  damping 
mechanism.  With  multiple  particles,  the  number  of  interparticle  and  particle-wall  impacts 
becomes  extremely  large  after  a  short  period  of  time.  Further,  it  is  difficult  to  describe  each 
impact  and  all  sources  of  mechanical  loss  accurately.  Thus,  there  are  fundamental  tradeoffs 
between  modeling  of  individual  particles  and  modeling  of  the  group  effects  of  the  collection  of 
particles.  There  is  also  a  need  to  correlate  models  with  experiments,  anchoring  the  analysis  to 
measured  data  from  typical  test  cases. 

The  remainder  of  this  report  introduces  basic  single-  and  multiple-article  damper  configurations 
and  considers  the  effects  of  variation  in  different  parameters  on  overall  performance.  Several 
systems  were  used  to  measure  particle  damping  performance,  and  these  are  described  next.  The 
following  sections  discuss  the  analytical  modeling  in  more  depth,  providing  theoretical 
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background.  Test  results  and  test-analysis  correlation  are  described  for  beam  and  chassis 
structures.  Measurements  for  the  proof  of  concept  structure  are  presented.  The  report  ends  with 
a  description  of  a  related  commercial  application  of  the  particle  damping  that  illustrates  its 
value  in  a  high-temperature  environment,  and  including  conclusions  and  recommendations  for 
further  research. 
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3  Background 

MPID  technology  has  been  derived  from  single-particle  impact  technology,  a  technology  that 
has  been  in  use  for  over  a  half  a  century.  Literature  typically  distinguishes  particle  damping 
from  impact  damping  based  on  the  number  and  sizes  of  the  auxiliary  masses  or  particles  in  a 
cavity.  As  shown  in  the  system  models  of  Figure  1,  MPID  is  a  term  that  is  used  to  imply 
multiple  auxiliary  particles  of  small  sizes  in  a  cavity,  while  impact  damping  is  a  term  that 
implies  a  single  auxiliary  mass  or  particle  in  a  cavity. 


Figure  1.  Single  Impact  Damper  and  Multiparticle  Impact  Damper  Schematics 

One  of  the  earliest  practical  applications  of  particle  damping  was  a  “bean  bag”  filled  with  lead 
shot  used  on  highway  road  signs.  In  recent  years,  numerous  analytical  and  experimental  studies 
have  been  performed  on  the  use  of  multiple  and  single  particle  impact  damping  concepts  for 
specific  applications.  A  commercial  sports  product  has  even  been  introduced  using  MPID  [1]. 

Impact  dampers,  single  or  multiple  particle,  exhibit  highly  non-linear  behavior,  introducing 
damping  to  a  system  through  several  mechanisms,  which  include  plastic  deformation,  external 
and  internal  friction,  and  momentum  transfer.  The  predominate  mechanism  by  which  a  damper 
dissipates  energy  from  a  vibrating  system  depends  on  certain  damper  parameters.  Previous 
studies  have  addressed  how  the  damping  mechanism  varies  with  cavity  fill  ratio,  coefficient  of 
restitution,  and  material  hardness  [2];  however,  previous  studies  have  failed  to  generate  a 
comprehensive  set  of  test  parameters  and  analytical  procedures  for  characterizing  the  response 
of  structural  systems  subject  to  MPID  under  dynamic  loads. 

MPID  provides  at  least  two  advantages  over  single-particle  impact  damping.  First,  the  peak 
level  of  impulsive  force  imparted  to  a  cavity  during  each  impact  is  reduced,  so  that  particle 
motion  is  less  likely  to  produce  undesirable  acoustic  noise  or  wear.  Second,  a  multiparticle 
impact  damper  is  easier  to  tailor  to  applications  in  which  spatial  constraints  for  mounting  a 
damper  are  severe,  i.e.,  damping  achieved  through  motion  of  multiple  particles  may  require 
much  smaller  cavity  volumes  than  those  required  in  achieving  damping  through  single  particle 
motion. 

The  most  extensive  body  of  research  published  on  MPID  concerns  damping  by  tungsten 
particles  introduced  to  a  low-frequency  single-degree-of-freedom  system  under  random 
excitation  [2].  This  work  showed  that  struetural  response  amplitudes  could  be  reduced  by 
nearly  50%,  depending  on  how  damper  parameters,  such  as  mass  ratio,  cavity  fill  ratio,  and 
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cavity  size,  were  construed.  Plots  in  Figure  2  show  how  particle  sizes  were  related  to  damping 
in  a  vibrating  system.  Amplitude  attenuation  is  displayed  relative  to  the  non-damped  system 
configuration,  with  the  larger  normalized  X-dimension  and  Y-dimension  representing  larger 
cavities. 


Figure  2.  Influence  of  Particle  Size  and  Cavity  Fill  Ratios  on  the  Effectiveness  of  MPID 

For  high  temperature  applications,  MPID,  unlike  viscoelastic,  vitreous  enamel,  or  friction 
damping,  offers  the  possibility  of  a  solution  that  is  rugged,  reliable,  and  simple  to  implement. 
Numerous  new  ceramic  materials  are  available  for  fabrication  of  particles,  as  shown  in  Figure 
3. 


Temperature,  °C 

Figure  3.  Ceramic  Material  Hardness  as  a  Function  of  Temperature 
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These  ceramics  and  others  offer  resistance  to  temperature,  corrosion,  or  thermal  aging.  To  date, 
however,  the  use  of  MPED  methods  is  rare:  without  a  rigorous  design  methodology,  application 
of  such  methods  has  been  a  time-intensive  and  costly  trial-and-error  process.  The  motivation 
for  developing  the  design  methodology  is  clear. 

Primary  parameters  have  been  identified  for  their  key  role  in  any  particle  damping  mechanism. 
Parameters  include:  damper  cavity  size,  particle  size,  shape,  material  properties,  and  damper 
cavity  material. 

Previous  research  concerning  particle  damping  considers  volume  fill  ratio,  the  ratio  of  cavity 
volume  occupied  by  particles  to  total  cavity  volume,  as  a  parameter  for  characterizing  damping 
capability.  CSA’s  work  takes  this  research  one  step  further,  recognizing  that  one  change  in 
volume  fill  ratio  can  be  achieved  several  ways — by  varying  particle  size,  number,  and/or  shape 
or  by  varying  cavity  size  or  shape.  Correspondingly,  one  change  in  volume  fill  ratio  does  not 
necessarily  result  in  one  unique  change  in  observed  damping  capability.  Thus,  fill  ratio  may  be 
characterized  in  terms  of  other,  more  basic,  parameters. 

Two  mechanisms,  dissipation  of  energy  through  impact  (heat  and  noise)  and  dissipation  of 
energy  through  friction,  have  been  identified  as  sources  of  particle  damping.  The  degree  to 
which  each  mechanism  participates  in  total  damping  capability  depends  on  the  primary 
parameters. 
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4  General  Characteristics  of  Particle  Damping 

As  a  precursor  to  developing  a  design  methodology,  CSA  and  UDRI  determined  that  the 
damping  mechanisms  of  particle  damping  needed  to  be  better  understood.  In  an  effort  to 
catalog  and  characterize  particle  damping,  many  experimental  tests  were  undertaken. 

Particle  damping  is  amplitude  dependent.  For  this  reason,  most  standard  testing  approaches  for 
determining  the  damping  of  a  system  cannot  be  applied  directly.  It  was  thus  necessary  to 
develop  test  approaches  that  could  properly  capture  the  particle  damping  behavior.  While  the 
developed  methods  are  useable  on  any  structure  to  which  it  is  desired  to  add  damping,  for 
comparison  purposes,  they  were  all  performed  on  a  simple  cantilevered  beam.  Figure  4  shows 
configurations  for  room-  and  elevated-temperature  testing.  Note  that  the  test  on  the  left  side  of 
Figure  4  was  configured  for  forced  excitation  and  the  elevated-temperature  test  on  the  right  was 
configured  for  transient  ring-down.  The  response  was  measured  by  an  accelerometer  at  or  near 
the  cavity. 


V  ■  ■ 


The  following  sections  will  compare  the  performance  of  single  particle  dampers  to 
multiparticle  impact  dampers. 

4. 1  Results  for  Single-Particle  Impactor  Configurations 

In  this  section,  representative  experimental  results  for  single  particle  impactor  configurations 
are  reviewed,  primarily  in  terms  of  transient  ring  down  event  testing.  In  Figure  6  the 
experimentally  measured  impact  initiated  transient  ring  down  of  the  baseline  (red)  and  a 
representative  single  particle  impactor  (blue)  are  compared  in  terms  of  the  raw  data  (left)  and 
the  postprocessed  damping  versus  amplitude  curves  (right).  The  triangular  look  of  the  raw  time 
data  of  the  damped  ring  down  is  typical  of  most  observed  single  particle  impactor  dominated 
events.  Also  evident  is  the  sharp  turn  off  of  the  impactor  effectiveness  around  1.75  seconds,  to 
this  region,  the  impactor  becomes  more  and  more  effective,  causing  a  rapid  increase  in 
effective  damping,  until  the  response  of  the  beam  is  reduced  enough  that  the  damper  essentially 
turns  off  and  the  structures  ring  down  transistions  to  its  baseline  behavior. 


Figure  6.  Time  and  Postprocessed  Damping  versus  Amplitude  Curves  for  a 
Representative  Single  Particle  Impact  Damper  Configuration 


The  postprocessed  results  on  the  right  of  Figure  6  are  generated  as  follows:  The  envelope  of  the 
ring  down  (extracted  using  the  Hilbert  transform)  is  taken  and  fitted  with  a  spline  fit  over 
multiple  peicewise  continous  discrete  regions  where  the  spline  fits  are  forced  to  have  the  same 
slope  as  their  neighboring  spline  fits  at  intersections.  The  spline  fits  can  then  be  used  to 
calculate  local  slope  which  itself  is  directly  related  to  the  damping  at  that  amplitude.  Due  to  the 
nature  of  the  implementation  of  the  fitting  code  (i.e.,  the  assumption  of  piece  wse  continuous 
nature  of  the  ring  down),  the  transition  point  for  impactor  is  ill-conditioi>iej2(ti  and  is  thus 
typically  not  shown. 

While  the  data  of  Figure  6  is  only  for  one  ring  down  case,  similar  data/trends  were  seen  in  most 
other  single  particle  impactor  ring  downs.  For  example.  Figure  7  shows  and  compares  two 
different  views  of  the  damping  versus  amplitude  trends  for  five  different  ring  downs  where  the 
only  variable  was  capsule  length.  While  all  relevant  data  sets  are  not  shown,  several  trends 
were  consistent  across  most  data  sets,  to  general  it  can  be  said  that 

•  the  longer  the  capsule,  the  higher  the  amplitude  at  which  the  damping  turned  off 
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the  longer  the  capsule  the  higher  the  peak  damping  (generally) 
the  more  the  added  mass,  the  higher  the  damping. 
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Figure  7.  Two  Additional  Views  of  a  Typical  Impactor  Ring  Down  Data  for  Various 
Capsule  Lengths  (the  Same  Data  Are  Shown  at  Slightly  Different  Angles) 


The  first  point  makes  intuitive  sense,  as  the  capsule  length  (for  a  fixed  particle  diameter)  sets 
the  rattle  space  and  once  the  response  level  is  lower  than  the  rattle  space,  the  particle  can  no 
longer  hit  both  walls  in  an  optimal  manner,  and  damping  rapidly  dies  down.  The  third 
conclusion  must  be  tempered  by  the  fact  that  if  damping  results  are  normalized  by  the  amount 
of  added  mass,  it  turns  out  that  the  smaller  particles  yielded  consistently  better  performance  in 
terms  of  damping  added  per  unit  added  mass.  Best  results  ranged  from  2.8  to  6.38  percent 
zeta/gram  for  these  configurations.  Best  results  were  generally  at  0.7  or  0.8  in  capsule  length. 


4.2  Results  for  Multiple  Particle  Configurations 

In  this  section,  typical  results  for  multiple-particle  damper  configurations  are  shown.  As  an 
example,  Figure  8  shows  the  representative  postprocessed  damping  versus  amplitude  for 
transient  ring  down  test  results  corresponding  to  three  different-sized  tungsten  carbide  particle 
sweep  test  sets.  The  data  are  plotted  versus  added  mass  in  order  to  give  a  fair  comparison 
across  test  sets.  With  the  exception  of  the  single-particle  impact  dominated  behavior  of  the  1/8- 
inch  diameter  spheres  at  low  fill  ratios,  damping  generally  increased  with  added  mass.  In  the 
1/16-inch  and  3/32-inch  test  data,  there  is  also  indication  that  after  a  certain  amount  of  added 
mass  (and  for  a  fixed  capsule  size  fill  ratio),  the  added  mass  is  not  adding  more  damping.  This 
makes  sense  for  the  cylindrical  capsule  tested.  Also  note  the  gradual  turn  off  behavior  at  low 
amplitudes  obvious  in  most  data,  especially  the  1/16-inch  and  3/32-inch  test  eases.  Note  also 
how  particle  size  can  play  a  role  in  damper  response  (i.e.,  the  3/32-inch  case  has  lower  damping 
levels  than  either  the  1/16-inch  or  1/8-inch  test  cases,  despite  having  equivalent  added  masses). 
In  this  case  the  effect  was  traced  to  packing  issues  relative  to  the  cavity  size.  Peak  damping 
levels  of  about  1.2  percent  zeta  were  achieved  in  the  best  cases.  This  is  significantly  lower  than 
the  levels  seen  for  the  single  particle  impactors. 
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Figure  8.  Representative  Results  for  Tungsten  Carbide  Sweeps  with  Different  Particle 
Diameters  (1/16-inch,  3/32-inch,  and  1/8-inch)  and  Side  Views  of  the  Same  Results 


Figure  9  depicts  a  representative  forced  dwell  sine  sweep  response,  for  a  MPID  attached  on  a 
cantilever  beam  that  is  excited  at  the  base  by  a  electromagnetic  shaker.  The  observed  transfer 
functions  illustrate  the  so-called  saddleback  phenomenon.  In  this  phenomenon,  as  excitation 
amplitude  is  increased,  the  particle  damped  system  becomes  more  effective  and  goes  through  a 
region  of  optimal  damping  (here  roughly  at  .25  V  drive  voltage)  and  then  proceeds  to  degrade 
in  performance.  Similar  amplitude  dependencies  are  seen  in  the  transient  ring  down  data. 
Finally,  Figure  10  depicts  typical  FFT  response  of  transient  ring  down  tests,  illustrating  how  the 
particle  damping  distributes  energy  from  the  damped  modes  into  other  spectral  bands.  The  level 
of  induced  spread  experimentally  depends  on  the  number  of  particles,  etc.  Clear  evidence  of 
tonal  or  harmonic  phenomena  is  obvious. 
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Eight  0.125  inch  tungsten  carbide  particles,  0.19  len  x  0.5  dia  cavity 


Figure  9.  Forced  Sine  Sweep  Test  Results  for  MPID  Illustrating  the  Saddleback 

Phenomena 


Figure  10.  Representative  Fast  Fourier  Transforms  (FFT)  Showing  Measured  Frequency 
Content  for  Three  Test  Cases  "  Frequency,  hz” 
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5  Effect  of  Variation  in  Damper  Parameters 

Primary  parameters  have  been  identified  for  their  key  role  in  any  particle  damping  mechanism. 
Parameters  include:  damper  cavity  size  and  shape,  particle  size,  particle  shape,  material 
properties,  and  number  of  particles.  Other  design  parameters,  such  as  the  volume  fill  ratio,  may 
be  described  as  a  function  of  the  primary  parameters. 

5.1  Particle  Size 

Figure  11  compares  the  amplitude-dependent  damping  for  two  experimental  data  sets  where 
added  mass  is  varied  with  all  test  parameters  (test  object,  capsule  shape/size/material,  particle 
shape,  added  mass  and  volume  for  each  comparable  test,  particle  material  type/density,  etc.) 
identical  with  the  exception  of  particle  size.  In  order  to  keep  the  added  mass  for  comparable 
test  constant,  a  proportionally  larger  number  of  particles  for  the  smaller  type  were  added  in 
order  to  keep  both  added  mass  and  volume  constant.  As  usual,  results  are  given  for  a  consistent 
scale.  The  results  for  the  smaller  particle  size  and  hence  greater  number  of  particles  are  shown 
on  the  left,  while  the  larger  particle  size  and  hence  fewer  actual  particles  are  shown  on  the  right. 
The  rough  edges  on  the  right  plot  are  due  to  the  fitting  routines  used  to  generate  the  surface  plot 
shown  and  do  not  represent  actual  data. 

Several  key  differences  are  noticeable.  The  most  important  is  that  the  response  amplitude  at 
which  peak  damping  occurs  can  be  shifted  through  proper  particle  size  selection.  For  this 
particular  cavity  type,  larger  particles  favor  lower  amplitude  (peaks  from  0  to  -5  dB)  while  the 
smaller  particles  favor  the  higher  amplitudes  (+3  to  +5  dB).  The  shape  of  the  peaks  are  also 
different  in  that  smaller  particles  yield  narrower  but  higher  peak  damping  values.  Finally,  the 
smaller  particles  are  seen  to  give  better  added  damping  at  higher  amplitudes  where  as  the  larger 
particles  yield  almost  none  (compare  the  20  to  30  dB  range).  These  trends  we  observed  across  a 
wide  range  of  material  type. 


X  10"*  X  10'' 


Response  Amplitude,  dB  re  1G  Response  Amplitude.  dB  re  1G 

Figure  11.  Comparison  of  Amplitude-Dependent  Damping  versus  Equivalent  Added  Mass 
for  Two  Different  Particle  Sizes  with  All  Else  Being  Constant 


12 


Total  particle  mass  appears  to  have  a  fairly  significant  effect  on  damping  for  both  single 
particle  and  multiple  particle  dampers.  Increasing  the  mass  tends  to  increase  damping. 
Conversely,  decreasing  the  mass  tends  to  reduce  damping.  As  might  be  expected,  changes  in 
the  total  particle  mass  can  lead  to  a  fairly  significant  shift  in  the  frequency  of  peak  response. 

5.2  Dependence  of  Observed  Behavior  on  Orientation  with  Gravity 

Figure  12  compares  the  variation  in  achieved  damping  for  nearly  identical  particle  damping 
treatments  in  a  fixed  capsule.  All  major  test  parameters  are  identical  (i.e.,  test  object,  capsule 
type/size,  material,  particle  size,  etc.)  except  for  the  orientation  of  the  test  object  relative  to  the 
local  quasi-static  acceleration  field.  As  indicated  by  the  small  cartoons  on  each  figure,  the  plot 
on  the  left  corresponds  to  the  case  where  the  excitation  is  not  aligned  with  the  local  quasi-static 
acceleration  field,  and  the  one  on  the  right  corresponds  to  the  aligned  test  condition.  There  are 
several  distinct  differences  in  behavior. 

For  example,  the  damping  for  the  aligned  direction  case  achieves  a  higher  added  damping  (2.44 
percent  zeta  versus  2.17  percent)  but  this  is  not  the  whole  story.  Perhaps  the  most  striking 
difference  is  how  consistently  and  sharply  the  damping  turns  off  once  the  response  levels  drop 
below  1  G  (or  0  dB  on  the  plots)  for  the  aligned  case.  Another  more  subtle  difference  is  for  the 
best  performing  damping  treatment  of  each  case  the  response  amplitude  range  over  which 
damping  is  added  at  the  2  percent  damping  is  about  identical  in  each  case.  At  1.5  percent  level, 
the  not  aligned  case  is  favored  slightly,  and  favored  quite  strongly  at  the  1  percent  damping 
level.  Also  the  extracted  ring  down  profiles  versus  added  mass  are  more  consistent  for  the 
aligned  case. 

A  key  conclusion  that  can  be  drawn  from  these  results  is  that  the  direction  in  which  the  local 
quasi-static  acceleration  field  acts  can  play  a  large  roll  in  determining  at  which  response 
amplitude  particle  damping  will  turn  off  and  it  also  helps  set  the  level  of  damping  that  can  be 
achieved. 


Not  Aligned  Aligned 


Figure  12.  Comparison  of  Amplitude-Dependent  Damping  versus  Added  Mass  for  Two 
Different  Orientations  of  the  Test  Object  Excitation  versus  Local  Quasi-Static 

Acceleration  Field 
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5.3  Dependence  on  Material  Parameters 

Figure  13  compares  the  amplitude-dependent  damping  for  two  experimental  data  sets  where  the 
number  of  added  particles,  and  hence  added  mass,  is  varied  with  all  test  parameters  (test  object, 
capsule  shape/size/material,  particle  size/shape,  number  of  particles  for  each  comparable  test, 
material  type/density,  etc.)  being  identical  with  the  exception  of  a  relatively  minor  variation  in 
tested  material  alloy.  The  material  on  the  left  was  SS-315.  The  material  on  the  right  was  SS- 
304.  While  they  both  show  similar  trends  for  increased  performance  with  added  particle 
count/mass,  the  one  on  the  left  is  clearly  superior  in  that  it  achieves  smoother  and  higher 
damping  levels  for  all  numbers  of  added  particle.  The  amplitude  of  peak  damping  also  does  not 
seem  to  vary  as  much.  The  high  amplitude  behavior  (i.e.,  above  20  dB)  for  both  types  is 
similar,  as  is  the  fact  that  they  turn  off  around  0  dB. 


Figure  13.  Comparison  of  Amplitude-Dependent  Damping  versus  Added  Mass  for  Two 
Experimental  Data  Sets  for  the  Same  Base  Material  but  Slightly  Different  Alloys 
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6  Test  System  Overview 

A  wide  range  of  test  rigs  (eight  basic  types)  were  developed  to  experimentally  investigate 
various  aspects  of  particle  damping  behavior.  The  name  and  general  purpose  of  each  of  the  test 
systems  is  listed  in  Table  1.  Pictures  of  each  of  the  test  system  are  gathered  in  Figure  14. 
Highlights  of  the  test  results  for  each  test  system  are  discussed  in  the  following  sections. 


Table  1.  Test  Systems 


Basic  Test  Rigs 

Name 

Purpose 

Basic  Beam  Tester 
Vertical  Config. 
Horizontal  Config. 

Measure  basic  particle  damping  performance  over  a  wide  range  of 

•  Particle  materials/shapes/sizes/fill  ratios 

•  Capsule  sizes/shapes/materials 

•  Orientation  of  excitation  versus  gravity 

•  Excitation  methods. 

“Big  Beam” 

Scaled  up  version  of  basic  beam  tester.  Used  to  experimentally  check 
scaling  laws  for  modal  mass  ratio,  particle  and  capsule  sizes,  etc.  Also 
used  to  study  the  effects  of  capsule  wall  thickness. 

Mass-Spring-Damper 
(MSD)  Approximation 

Study  mass  ratio,  stiffness,  and  inherent  damping  effects  in  a 
controlled  manner.  Also  ensures  linear  motion. 

Temperature  Sensitivity 

Temperature  Chamber 
Sweeps 

Measure  the  effect  of  range  of  intermediate  temperatures  (-94  to  +365 
Tl  on  the  behavior  of  particle  damping 

Kiln  Test  Rig 

Measure  the  effect  of  high  temperatures  (70  to  850  T)  on  the  behavior 
of  particle  damping. 

Higher  Order  Structure  Test  Rigs 

Panel  Test  Rig 

Demonstrate  particle  damping  can  work  on  more  complicated 
structures  and  study  design  approaches. 

High-Temp  Material 
Test  Rig 

Demonstrate  particles  damping  can  work  on  complicated  structures 
and  at  high  temperature. 

Other 

Duration  Tester 

Begin  preliminary  measurements  of  durability  of  particle  damping 
treatments. 
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Figure  14.  Overview  of  the  Various  Test  Systems,  Left  to  Right  Working  Down  from 
Upper  Left:  Basic  Beam  (Horizontal  Configuration),  Oven,  Panel,  Kiln,  Big  Beam, 
Duration  Test  Rig,  Mass-Spring-Damper 


6. 1  Basic  Beam  Tester 

Most  of  the  experimental  studies  of  particle  damping  behavior  were  performed  with  the  test 
object  shown  in  Figure  15.  It  consisted  of  a  well-clamped  cantilever  beam.  Per  the  standard  test 
approach  experimental  displacement  or  acceleration  data  is  recorded  over  the  ring  down  time. 
A  Hilbert  transform  is  used  to  capture  the  envelope  of  this  data.  The  envelope  is  subsequently 
fit  with  splines  over  multiple  discrete,  but  piecewise  continuous,  regions.  These  spline  fits  can 
then  be  used  to  calculate  the  local  slope  of  the  ring  down,  which  is  related  to  the  damping  at 
that  amplitude.  This  approach  is  similar  to  the  log  decrement  method,  but  is  less  sensitive  to 
noise. 

Well  over  2,000  individual  ring  down  tests  of  different  combinations  of  particles,  capsules,  and 
orientation  relative  to  gravity  have  been  performed  with  this  test  setup.  Representative  results 
from  some  of  these  tests  are  shown  in  Figure  16.  The  data  compares  the  differences  in  behavior 
for  the  different  amounts  of  particles  in  the  same  capsule.  The  test  setup  has  been  used  to  study 
the  effects  of  the  following. 

Particles 

•  Different  materials 

•  Different  heat  treats  on  otherwise  similar  materials 

•  Particle  size 

•  Particle  shape 

•  Particle  fill  ratio 

•  Number  of  particles  (single  versus  multiple) 

Capsules 

•  Capsule  diameter 

•  Capsule  length 

•  Capsule  material  (limited) 

Orientation  of  excitation  relative  to  local  gravity 

•  Aligned 

•  Not  aligned 

Excitation  approaches 

•  Air  pulses  (insufficient  excitation) 

•  Electrodynamic  shaker  (used  for  sine  sweeps  and  forced  response  only) 

•  Manual  displacements  or  tapping. 
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Figure  15.  Experimental  Test  Setup  Used  for  Basic  Particle  Damping  Characterization  (in 
the  Excitation  Aligned  with  Gravity  Orientation) 


Figure  16.  Experimental  Results  Comparing  Amplitude-Dependent  Damping  for  Multiple 
Particle  Dampers  with  (left)  Numerous,  Smaller  Particles  versus  (right)  Fewer,  Larger 

Particles 


6.2  Big  Beam 

The  objective  of  this  test  hardware  was  to  gather  data  with  much  larger  dampers  to  determine 
if/how  the  data  for  smaller  capsules  scaled  to  that  of  the  larger  particle  dampers,  and  if  larger 
dampers  could  produce  similar  performance  as  that  of  the  smaller  ones.  The  experimental 
apparatus  consisted  of  the  beam,  two  identical  particle  damper  assemblies,  and  the 
instrumentation  and  data  acquisition  system.  The  host  structure  used  was  a  large,  vertically 
oriented,  cantilevered  beam.  The  beams  dimensions  were  4  by  1  inches  with  a  free  length  of 
42-9/16  inches.  Two  independent  capsule  holders  could  accommodate  capsules  that  ranged 
from  1.1, 1.6,  and  2.1  inches  in  inner  diameter,  and  lengths  from  0.6  to  4.1  inches  long. 

An  initial  sweep  with  the  largest  ID  capsule  and  five  different  particle  diameters  revealed  that 
the  longest  capsule  and  the  second  largest  particle  size  yielded  the  best  results.  More  detailed 
tests  with  this  capsule  revealed  that,  unlike  the  smaller  damper  capsules  which  generally  had 


peak  damping  about  50  percent  full,  these  capsules  had  best  damping  when  nearly  100  percent 
full.  Also,  after  crossing  the  50  percent  fill  ratio,  the  behavior  of  the  damping  changed  slightly, 
as  shown  in  Figure  17.  The  plot  in  Figure  17  shows  a  specific  case  where  the  number  of 
particles  was  increased  from  50  to  176  particles  using  7/16-inch  stainless  steel  balls.  The  cavity 
was  essentially  100%  full  with  176  balls. 

The  larger  physical  size  of  the  damper  capsules  also  allowed  the  initial  study  of  the  effects  of 
different  wall  and  cap  thickness.  For  the  limited  tests  performed,  it  was  seen  that  both  wall 
thickness  and  cap  thickness  can  influence  the  observed  damping,  actually  improving  it  in  the 
measured  test  cases  (6.5%  peak  damping  for  thick  case,  9.0%  peak  for  the  thin/thin  case). 


e/ia^01  Large  Particle  Damper  Results 
4,r  X  2.1  Thick  Wall  Capsule 
y/ic  SS  Bans 


Ringrtown  Amplitude,  dB  re  1  G 

Figure  17.  Representative  Test  Results  for  the  Big  Beam  Tester  Where  Increase  in  the 
Number  of  Particles  (Up  to  Nearly  100  Percent  Full)  Showed  Increase  in  Damping 

6.3  Mass  Spring  Damper  Approximations 

Using  modified  hardware  taken  from  previous  projects  and  additional  custom  test  components, 
a  test  apparatus  which  served  as  an  engineering  approximation  of  the  standard  linear,  single- 
degree-of-freedom  (SDOF)  MSD  system  was  created.  Using  this  test  apparatus  and  a  fixed 
particle  damping  treatment,  it  was  then  possible  to  individually  study  the  effects  of  variations  in 
stiffness,  mass,  and  mass  ratios  as  well  as,  to  a  limited  extent,  see  how  particle  damping 
performed  in  conjunction  with  other  damping  treatments. 

This  test  setup  was  useful,  as  the  majority  of  the  experimental  tests  described  in  this  document 
were  on  cantilevered  beams  or  more  complex  structures.  This  made  it  difficult  to  conclusively 
decouple  the  effects  of  variations  in  mass  ratio,  stiffness,  etc.  There  was  also  the  slight  added 
complication  of  the  fact  that,  at  least  for  the  cantilever  beam  cases,  the  path  of  the  capsule 
motion  is  not  truly  linear  but  follows  a  curved  profile  set  by  the  length  of  the  beam.  This  adds 
another  discrepancy  between  experimental  work  and  explicit  modeling  results.  Having  a  linear 
spring  element  reduces  this  discrepancy  and  provides  the  ability  to  vary  mass  and  stiffness 
values. 

The  MSD  test  hardware  is  shown  in  Figure  18.  It  consisted  of  a 
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Standard  test  capsule  holder  and  spiral  flexure  end  clamp 

Two  or  more  spiral  flexures 

Spacer  housing  and  spiral  flexure  end  caps 

Central  spacing  element 

Spiral  flexure  end  clamp 

Accelerometer  mount 

Accelerometer 

Additional  moving  mass  elements 


Figure  18.  Mass-Spring-Damper  Approximation  Setup 


Stiffness  Elements 

Stiffness  was  provided  by  spiral  flexure.  Each  flexure  provides  roughly  12,400  N/m  of  stiffness 
(equivalent  to  70.8  Ibf/in).  The  total  system  level  stiffness  is  then  N  x  12,400  N/m  where  N  is 
the  number  of  implemented  flexures.  The  minimum  number  for  N  is  2,  maximum  is  10  but  6 
was  the  maximum  which  yielded  good  test  results 

Moving  Mass 

The  basic  structure  of  the  MSD  hardware  created  an  inherent  baseline  moving  mass  of 
approximately  170  grams.  In  order  to  perform  tests  where  the  mass  and  stiffness  increased  in 
equal  amounts  (hence  keeping  the  natural  frequency  the  same  but  changing  the  mass  ratio)  it 
was  necessary  to  have  custom  mass  elements  made.  They  bolted  onto  the  outside  of  the  MSD 
for  ease  of  integration.  Four  additional  masses  corresponding  to  the  4  additional  stiffness  states 
(N  =  3, 4,  5, 6)  were  made.  The  mass  increment  turned  out  to  be  82.8  grams. 

Inherent  Damping 

Aerodynamic  damping  was  one  source  of  damping  for  the  two-flexure  case.  When  more  than 
two  flexures  were  used,  interlaminar  friction  between  neighboring  flexures  contributed 
significantly  more  damping.  Inherent  damping  levels  of  approximately  0.4  percent  to  2.0 
percent  zeta  were  experimentally  observed.  At  first,  these  large  levels  of  inherent  damping 
were  feared  to  create  too  high  of  a  ‘damping  noise  floor’  but  later  it  was  seen  that  this  was  not 
the  case. 


Particle  Damping  Treatment 


In  order  to  enable  fair  comparisons  across  all  test  configurations,  a  single  particle  damping 
treatment  was  used.  It  consisted  of  7  by  9/64  inch  diameter  bores  in  a  titanium  test  capsule  with 
two  caps.  Free  internal  length  was  5/32-inch. 

Due  to  time  constraints,  the  test  matrix  shown  in  Table  2  was  only  implemented  twice.  Once 
with  the  already  described  damping  treatment,  and  then  once  again  with  the  particles  removed. 
Over  50  individual  tests  were  performed. 

Table  2.  Test  Matrix 


Added  Mass 

Flexures 

2 

3 

4 

5 

6 

# 

Stiffness 

169.3g 

243.9  g 

329.7  g 

415.5  g 

496.1  g 

2 

24,800  N/m 

1.00 

0.82 

0.71 

0.63 

0.58 

3 

37,200  N/m 

1.22 

1.00 

0.87 

0.77 

0.71 

4 

49,600  N/m 

1.41 

1.15 

1.00 

0.89 

0.82 

5 

62,000  N/m 

1.58 

1.29 

1.12 

1.00 

0.91 

6 

74,400  N/m 

1.73 

1.41 

1.22 

1.10 

1.00 

An  enforced  displacement  was  used  upon  release  to  impart  a  starting  acceleration  to  the  test 
object.  The  resulting  free  decay  was  then  measured  using  an  accelerometer  mounted  collinearly 
with  the  mass. 

Peak  damping  is  normally  expected  to  occur  when  the  particles  hit  both  ends  of  the  capsule  in 
phase,  2  times  per  cycle.  This  occurs  at  a  fixed  length/displacement  which  implies  that  as  the 
resonant  frequency  changes,  the  acceleration  level  at  which  peak  damping  occurs  should  be 
changing.  This  was  not  observed.  In  fact,  the  peak  damping  generally  occurred  around  the  same 
acceleration  value. 

A  comparison  of  the  trends  in  levels  of  intentionally  added  damping  versus  mass  ratio  for  the 
various  tested  stiffness  cases  (where  k  is  the  number  of  flexures)  is  given  in  Figure  19.  In  this 
plot  the  peak  observed  damping  minus  the  mean  of  the  observed  baseline  damping  are  plotted. 
As  expected  there  is  a  clear  decrease  in  added  damping  as  the  mass  ratio  increases.  As  there  is 
some  noise  in  the  data  it  is  hard  to  say  what  trend  the  decrease  specifically  follows.  The 
originally  hypothesized  strict  ratio  rule  is  indicated  with  the  dashed  lines  for  different  level  of 
initial  damping. 

Interesting  results/conclusions  from  these  tests  include  the  following. 

•  Damping  was  seen  to  experimentally  depend  on  the  ratio  of  the  damping  treatment  mass 
to  the  total  moving  mass. 

•  Amplitude  of  peak  damping  did  not  appear  to  vary  with  frequency. 
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Figure  19.  Representative  MSD  Test  Results 


6.4  Temperature  Chamber  Sweeps 

The  purpose  of  this  test  setup  was  to  enable  testing  of  variable  capsule  geometries  and  particle 
fills  in  the  standard  cantilever  beam  over  a  wide  temperature  range  (-70  to  +185  ”C,  -94  to  365 
T).  This  corresponds  to  a  255  °C  (459  T)  total  temperature  range.  It  was  intended  to 
complement  the  high-temperature  kiln  (described  in  Section  6.5)  in  that  it  enabled  testing  to 
very  cold  levels  well  below  room  temperature  and  tighter  temperature  control  at  elevated 
temperature  up  to  the  environmental  chamber’s  upper  limit. 

The  test  setup  consisted  primarily  of  the  standard  cantilever  beam.  As  the  test  beam  was  now 
inside  of  an  environmental  chamber,  an  electromechanical  push-pull  solenoid  system  was  used 
to  create  the  required  initial  displacement  and  release.  This  turned  out  to  provide  a  very 
repeatable  excitation. 

For  the  majority  of  the  performed  tests,  a  titanium  capsule  with  seven  through  holes  containing 
seven  tungsten  carbide  spheres  1/8-inch  in  diameter  was  used.  The  standard  data  acquisition 
system  and  postprocessing  approach  using  Hilbert  transforms  of  the  ring  down  was  applied. 
Representative  results  are  shown  in  Figure  20  where  damping  (percent  zeta)  is  plotted  versus 
amplitude  (dB  relative  to  1  G).  Comparative  baseline  curves  that  remain  below  0.4  zeta 
represent  tests  at  each  temperature  where  no  additional  damping  was  intentionally  added.  They 
are  shown  in  solid  colors.  The  particle  damping  measurements  for  the  given  temperatures  are 
dashed. 

While  there  is  some  variability  in  the  baseline  measurements  and  in  the  damped  measurements, 
the  seen  changes  do  not  correlate  to  changes  in  temperature,  and  are  believed  to  represent 
inherent  noise  in  the  test/data  extraction  process.  The  overall  scale  of  observed  damping  for  the 
intentionally  damped  cases,  versus  the  baseline  cases  (i.e.,  3  percent  intentional  versus  0.35 
percent  baseline  peak  damping)  is  significant.  The  key  conclusion  is  that  very  little  sensitivity 
to  temperature  was  seen  over  the  test  range. 


22 


Figure  20.  Representative  Baseline  and  Damped  Experimental  Test  Results  at  Various 
Temperatures  and  the  Test  Setup  Installed  in  the  Environmental  Chamber 


6.5  Kiln  Test  Rig 

The  ultimate  goal  of  this  testing  effort  was  to  see  if  the  damping  characteristics  of  an  impact  or 
particle  damper  would  change  with  temperatures  of  up  to  1,000  °F.  The  temperature  limitations 
of  transducers  prevented  the  experimenters  from  collecting  data  up  to  that  temperature. 
Reliable  data  was  collected  from  room  temperature  up  to  850  T. 

The  tube-furnace  shown  in  Figure  21  was  modified  to  house  a  stainless-steel  fixture  and 
cantilever  beam.  The  excitation  was  performed  using  a  spring-loaded  plunger  system. 
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Figure  21.  High-Temperature  Test  Rig 


Thermal  expansion  contributed  significantly  to  problems  and  irregularities  encountered  during 
initial  testing.  The  original  design  of  the  test  apparatus  had  to  be  modified  slightly  (Bellville 
washers)  in  order  to  limit  the  effects  of  thermal  expansion  of  the  hardware  on  the  measured 
data.  This  is  shown  in  Figure  22.  It  is  believed  that  most  of  these  problems  were  overcome. 
Once  these  issues  were  removed,  it  was  experimentally  demonstrated  that  particle  damping  can 
be  applied  to  applications  up  to  at  least  850  T.  In  some  cases,  some  sensitivity  to  temperature 
was  seen  (damping  actually  improved).  Representative  data  is  shown  below  for  a  multiple 
particle  cavity  setup  where,  again,  the  post  processing  approach  using  Hilbert  transforms  of  the 
ring  down  was  applied.  Representative  results  at  various  temperatures  are  shown  in  Figure  23 
where  damping  (percent  zeta)  is  plotted  versus  amplitude  (dB  relative  to  1  G).  While  the  data 
are  noisy,  there  are  no  clear  temperature  dependencies.  The  baseline  curve  is  a  room 
temperature  test  with  no  particles  in  the  cavity. 
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6.6  Duration  Tester 

The  objective  of  this  test  hardware,  shown  below  in  Figure  24,  was  to  determine  if  MPE)s 
undergo  any  significant  wear  and/or  changes  in  damping  characteristics  as  they  are  subjected  to 
many  cycles.  The  experimental  apparatus  consisted  of  the  particle  damper  itself,  a  cantilevered 
beam  host  structure,  excitation  device,  and  data  acquisition  system.  Two  different  particle 
dampers  were  tested.  The  testing  performed  revealed  that  the  dampers  did  undergo  significant 
wear  over  about  6  days  of  continuous  operation.  From  a  product  standpoint,  this  could  be  a 
serious  concern  and  design  challenge.  Fortunately,  the  wear  that  did  occur  tended  to  increase 
the  performance  of  the  damper. 


7  Theoretical  Background 

The  particle  damper  is  conceptually  a  relatively  simple  device.  However,  the  behavior  of  the 
particle  damper  is  highly  nonlinear  and  energy  dissipation,  or  damping,  is  derived  from  a 
combination  of  loss  mechanisms.  The  relative  effectiveness  of  these  mechanisms  changes 
based  on  various  system  parameters.  In  the  first  subsection,  loss  mechanisms  present  in  the 
particle  damper  are  reviewed. 

In  the  next  subsection,  the  incorporation  of  these  loss  mechanisms  into  the  mathematical  model 
is  discussed.  The  mathematical  model  has  been  developed  based  on  the  particle  dynamics 
method.  The  utility  of  the  particle  dynamics  method  is  derived  from  the  ability  to  simulate 
contact  interaction  forces  with  a  small  number  of  parameters  which  capture  the  most  important 
contact  properties.  One  of  the  critical  aspects  for  developing  an  accurate  mathematical  model  is 
the  selection  of  appropriate  force-displacement  relations. 

Lastly,  the  incorporation  of  the  particle  dynamics  method  into  an  explicit  finite  element  code  is 
discussed.  Contact  detection  algorithms  are  reviewed,  along  with  the  implementation  of  the 
force-displacement  relations  used  to  resolve  the  contact.  Sample  results  and  correlation  of 
predicted  results  with  measured  experimental  data  are  considered  in  Sections  4  and  5. 

7. 1  Particle  Damper  Loss  Mechanisms 

Passive  damping  techniques  attenuate  the  response  of  a  vibrating  structure  by  removing  a 
portion  of  the  vibratory  energy.  One  method  of  removing  the  vibratory  energy  is  to  transfer  the 
energy  to  a  secondary  system.  Dynamic  vibration  absorbers  function  in  this  manner. 
Oscillations  of  the  dynamic  vibration  absorber  are  created  with  a  fraction  of  the  energy  creating 
the  oscillations  of  the  primary  system.  This  mechanism  is  present  in  the  particle  damper  in  the 
form  of  the  momentum  transfer  which  occurs  when  the  particles  impact  the  cavity  ends. 

A  second  method  of  removing  the  vibratory  energy  is  to  dissipate  the  mechanical  energy  in  the 
form  of  heat  or  noise.  Viscoelastic  and  viscous  damping  dissipate  energy  in  the  form  of  heat 
created  when  the  viscoelastic  material  or  viscous  fluid  undergoes  shear.  Friction  damping  also 
dissipates  energy  in  the  form  of  heat  which  is  created  due  to  the  relative  motion  between  two 
contacting  surfaces.  In  particle  dampers,  friction  exists  when  relative  sliding  (or,  to  a  much 
lesser  degree,  rolling)  motion  occurs  between  individual  particles  in  contact  or  between 
particles  in  contact  with  the  cavity  walls.  Viscoelastic  behavior  also  is  present  in  the  particle 
dampers  due  to  impacts  with  nonunity  coefficients  of  restitution  (i.e.,  impacts  which  are  not 
perfectly  elastic). 

For  convenience,  the  loss  mechanisms  present  in  the  particle  damper  can  be  grouped  into 
“external”  and  “internal”  mechanisms.  “External”  mechanisms  involve  friction  and  impact 
interactions  between  the  particles  and  the  cavity.  “Internal”  mechanisms  involve  friction  and 
impact  interactions  between  the  individual  particles.  The  relative  effectiveness  of  these 
mechanisms  changes  based  on  various  system  parameters.  For  example,  at  low  excitation  levels 
or  for  relatively  long  eavity  lengths,  insufficient  energy  may  exist  to  create  impacts  between  the 
particles  and  the  cavity  ends.  Under  such  circumstances,  any  dissipation  would  be  due  solely  to 
any  sliding  or  rolling  friction  between  the  particles  and  the  cavity  or  interactions  between  the 
individual  particles.  For  higher  excitations  or  shorter  cavity  lengths,  however,  losses  due  to 
impacts  between  the  particles  and  the  cavity  ends  may  dominate  the  overall  dissipation. 
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Incorporation  of  the  various  loss  mechanisms  into  the  mathematical  model  will  be  performed 
through  the  definition  of  appropriate  force-displacement  relations. 

7.2  Force-Displacement  Relations 

In  the  following  paragraphs,  the  force-displacement  relations  used  to  account  for  the  contact 
interactions  between  the  individual  particles  and  between  the  particles  and  the  cavity  walls  are 
considered.  It  has  been  assumed  that  the  particle  dampers  consist  of  spherical  particles  of  the 
identical  material.  Consider  a  typical  impact  of  two  spherical  particles,  i  and  j,  with  radii  Ri  and 
Rj,  with  the  particle  centers  separated  by  a  distance,  dij,  as  shown  in  Figure  25. 


Figure  25.  Typical  Particle-Particle  Impact  Parameters 

These  two  particles  interact  if  their  approach,  a,  is  positive.  The  approach  can  be  defined 
according  to  the  following  equation 


o  =  (R,+Rj)-d„  (1) 

In  this  case,  the  colliding  spheres  feel  a  force,  F,  equal  to  the  following 

F  =  F^-n^+F^-n^  (2) 

where  and  n®  are  the  unit  vectors  in  the  normal  and  shear  directions,  respectively,  and  f'^ 
and  F®  refer  to  the  normal  and  shear  forces.  The  normal  force  acting  between  colliding  spheres 
can  be  broken  down  into  an  elastic  portion  and  a  dissipative  portion  as  follows: 

F"=F5)  +  FSs)  (3) 

The  following  paragraphs  discuss  the  elastic  and  dissipative  portions  of  the  normal  force,  the 
shear  force,  and  the  forces  due  to  particle-cavity  impacts. 
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7.2.1  Elastic  Portion  of  Normal  Force 

Expressions  for  the  elastic  portion  of  the  normal  force  are  based  on  Hertz’s  theory  of  elastic 
contact  [3-4].  Hertz  derived  expressions  for  the  stress  between  contacting  bodies  based  on  the 
basic  equations  of  elasticity  and  on  the  geometries  of  the  contacting  bodies.  More  specifically, 
Hertz’s  expressions  are  an  extension  of  the  work  of  Boussinesq  [5],  who  used  the  theory  of 
potential  to  establish  a  method  to  determine  the  stresses  and  displacements  resulting  due  to 
surface  tractions  on  an  elastic  half-space.  The  derivation  of  Hertz’s  theory  is  discussed  briefly 
in  the  following  paragraphs.  Complete  derivation  of  Hertz’s  theory  can  be  found  to  various 
levels  of  detail  in  a  number  of  classical  elasticity  texts  [6-8]. 

For  elastic  bodies  in  contact,  Hertz  determined  that  the  contact  area  is  generally  elliptical  and 
assumed  that  each  body  could  be  regarded  as  an  elastic  half-space  loaded  over  a  small  elliptical 
region  of  its  plane  surface.  This  assumption  implies  that  the  dimensions  of  the  contact  area 
must  be  small  compared  with  the  dimensions  of  each  body  and  with  the  relative  radii  of 
curvature  of  the  surfaces.  These  assumptions  are  necessary  to  ensure  that:  (1)  the  boundaries  of 
the  bodies  do  not  significantly  affect  the  stresses  in  the  contact  zone;  (2)  the  portions  of  the 
bodies  outside  of  the  contact  zone  can  be  roughly  approximated  as  an  elastic  half-space;  and, 
(3)  the  strains  remain  sufficiently  small  that  the  linear  theory  of  elasticity  applies.  Hertz  also 
assumed  that  the  contact  is  frictionless  such  that  only  normal  pressures  are  transmitted  through 
the  contact. 

For  two  elastic  bodies  in  contact.  Hertz  showed  that  the  total  normal  force  is  related  to  the 
approach  via  the  power  law 

F(^)=Ca-  (4) 

where  C  is  a  constant  which  depends  on  the  elastic  properties  of  the  materials  and  on  the  local 
curvature  of  the  bodies.  For  the  case  of  two  contacting  spheres  with  identical  properties,  a 
circular  contact  area  with  radius,  a,  results.  Hertz  showed  that  Equation  4  becomes 


where 


(R,R,) 

(Ri+Rj) 


(6) 


and  E  and  v>  are  the  Young’s  modulus  and  Poisson’s  ratio  of  the  spheres,  respectively.  The 
approach  and  contact  circle  radius  are  related  as  follows: 


(7) 


Equation  5  also  holds  for  two  impacting  spheres,  provided  that  the  duration  of  the  collision  is 
long  compared  with  the  first  fundamental  mode  of  vibration  in  the  spheres  [9]. 


7.2.2  Particle-Cavity  Relations 

In  addition  to  particle-particle  impacts,  particle-cavity  impacts  also  can  occur.  Saluena,  et  al. 
[10]  resolve  particle-cavity  impacts  by  using  a  fixed  layer  of  particles  along  the  cavity  walls.  A 
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similar  procedure  has  been  used  in  some  experimental  efforts  [11].  This  assumption  may  be 
appropriate  for  simulating  behavior  such  as  granular  flow  down  an  inclined  chute,  where  there 
may  exist  a  thin  layer  of  particles  which  essentially  remain  in  contact  with  the  cavity  wall. 
However,  particles  in  a  damper  often  may  collide  with  the  cavity  walls  rather  than  with  other 
particles.  As  a  result,  particle-cavity  force-displacement  relations  are  required.  Although  the 
complete  derivations  are  not  shown  here,  the  particle-cavity  force-displacement  relations  can  be 
formulated  by  modifying  the  particle-particle  force-displacement  relations  to  account  for  the 
material  properties  of  the  cavity  and  the  local  curvature.  For  simplicity,  it  has  been  assumed 
that  the  cavity  walls  are  flat  and  rigid.  Particle-cavity  relations  for  flat  cavity  walls  are  derived 
from  the  particle-particle  relations  by  assuming  the  radius  of  curvature  for  the  cavity  walls  goes 
to  infinity.  Note  that  this  assumption  is  valid  since  the  flat  surface  more  closely  represents  the 
Boussinesq  solution  for  the  surface  tractions  on  an  elastic  half-space  on  which  Hertz’s  theory  is 
based.  Under  the  assumption  of  rigid  walls,  the  effective  modulus  simply  becomes  the  particle 
modulus. 

7.2.3  Dissipative  Portion  of  Normal  Force 

For  damping  studies,  it  is  important  to  model  the  dissipative  portion  of  the  normal  force.  If  the 
granular  material  in  the  particle  dampers  were  to  behave  in  a  gas-like  manner,  it  would  be 
possible  to  resolve  particle-particle  and  particle-cavity  impacts  instantaneously  using  a 
coefficient  of  restitution  or  similar  parameter.  However,  the  granular  material  in  particle 
dampers  tends  to  behave  more  like  a  liquid  or  solid  in  that  enduring  contacts  occur.  As  a  result, 
the  use  of  a  coefficient  of  restitution  is  no  longer  appropriate  and  any  dissipation  must  be 
accounted  for  in  the  force-displacement  relations. 

Studies  on  the  coefficient  of  restitution  for  the  impact  of  two  spheres  have  demonstrated  that 
energy  is  dissipated  due  to  the  viscoelastic  behavior  of  the  sphere  material  [12].  Finite  element 
analyses  of  a  sphere  impacting  a  rigid  plate  have  been  performed  to  investigate  the  significance 
of  viscoelasticity.  An  axisymmetric  model,  shown  in  Figure  26,  has  been  used  for  the  analyses. 
The  model  contains  16,000  4-node  bilinear  axisymmetric  solid  elements  (CAX4  elements  in 
ABAQUS)  to  accurately  capture  the  contact  behavior.  Figure  27  shows  the  region  where  the 
spherical  particle  initially  contacts  the  rigid  wall. 


30 


Figure  26.  Axisymmetric  Model  of  Spherical  Particle  and  Rigid  Wall 


Figure  27.  Region  Where  Spherical  Particle  Initially  Contacts  Rigid  Wall 

To  validate  the  finite  element  model,  initial  analyses  have  investigated  the  impact  of  an  elastic 
sphere  of  radius  0.0625  inch  with  a  rigid  wall  and  have  been  compared  to  predicted  results 
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using  Hertz’s  theory.  The  sphere  is  given  steel  properties  with  an  elastic  modulus  of  30x10^ 
lbf/in^  a  density  of  0.283  lbf/in^  and  a  Poisson’s  ratio  of  0.33.  The  sphere  is  given  an  initial 
velocity  of  62.832  in/sec.  Figure  28  shows  predicted  load  and  displacement  curves  as  a  function 
of  time  from  both  the  finite  element  analyses  and  as  calculated  using  Hertz’s  relation.  As  seen 
in  the  figure,  the  hand  calculations  and  predicted  results  from  ABAQUS  compare  well. 


Figure  28.  Calculated  and  Predicted  Elastic  Behavior 

Viscoelastic  material  behavior  is  characterized  by  delayed  elasticity  (relaxation)  and  creep. 
Such  behavior  is  seen  in  many  materials  and  results  in  time-dependent  constituitive  equations. 
Linear  viscoelasticity  is  assumed,  and  the  strains  must  remain  small  and  the  principal  of 
superposition  must  apply. 

There  are  numerous  ways  to  model  a  linear  viscoelastic  material.  For  this  effort,  the 
viscoelastic  material  is  modeled  as  a  generalized  Maxwell  model  [13],  as  shown  in  Figure  29. 
This  model  includes  an  elastic  solid  (as  indicated  by  the  spring  element)  in  parallel  with  a 
number  of  Maxwell  fluid  elements  (as  indicated  by  elements  with  a  spring  and  dashpot  in 
series)  and  can  be  made  to  fit  relaxation  data  quite  accurately.  The  stress-strain  relationship  can 
be  expressed  using  a  relaxation  function,  ^(t),  which  specifies  the  stress  response  to  a  unit  step 
change  of  strain.  Similarly,  a  creep  compliance,  0(t),  which  specifies  the  strain  response  to  a 
unit  step  change  of  stress,  can  be  used.  Assuming  Poisson’s  ratio  remains  constant  in  time,  the 
relaxation  function  for  a  viscoelastic  material  modeled  as  shown  in  Figure  29  can  be  expressed 
as  follows: 


32 


(8) 


'F(t)=Eo  +Eie'‘'''  +E2e”'^'^  ...  +  E„e''''" 

In  Equation  8,  the  instantaneous  modulus  is  represented  as  the  sum  of  the  Ei,  and  the  elastic 
modulus  at  long  times  is  represented  by  Eq.  The  moduli,  Ei,  and  corresponding  time  constants, 
Ti,  can  be  determined  from  complex  modulus  material  properties,  such  as  those  obtained  by 
testing  in  accordance  with  ASTM  E756-93,  “Standard  Test  Method  for  Measuring  Vibration- 
Damping  Properties  of  Materials”  [14]. 


Figure  29.  Maxwell  Model  for  Viscoelastic  Behavior 

The  time  constants  in  Equation  8  are  used  to  account  for  the  behavior  of  the  material  during 
loading  over  various  time  periods  or,  equivalently,  at  various  frequencies.  From  Hertz’s  theory, 
the  time  of  maximum  compression  for  an  elastic  impact  of  two  spheres  can  be  calculated.  As  a 
first  approximation,  the  time  to  maximum  compression  can  be  thought  of  as  a  quarter  cycle  of 
sinusoidal  loading.  Thus,  approximate  loading  frequencies  can  be  calculated  based  on  the  time 
to  maximum  compression.  The  compression  times  for  the  impact  of  two  identical  elastic  steel 
spheres  at  various  impact  velocities  which  might  be  expected  in  a  typical  particle  damper  have 
been  calculated.  Results  from  these  calculations,  along  with  the  approximate  period  of  a 
loading  cycle  and  the  approximate  loading  frequency  are  given  in  Table  3. 


Table  3.  Approximate  Loading  Frequencies  for  Various  Impact  Velocities 


1  Impact  Velocity 
(in/s) 

Time  of  Maximum 

Compression  (s) 

Approximate  Period 

of  Loading  (s) 

Approximate  Loading 

Frequency  (Hz) 

6.28x10'^ 

1.72x10'^ 

6.90x10° 

1.45x10'' 

6.28x10'^ 

1.09x10'^ 

4.35x10° 

2.30x10'' 

6.28x10° 

6.86x10° 

2.75x10'° 

3.64x10'' 

6.28x10' 

4.33x10° 

i 

1.73x10'° 

5.77x10'' 
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6.28x10^ 

2.73x10'® 

1.09x10'^ 

9.15x10^ 

6.28x10^ 

1.72x10'® 

6.90x10'® 

1.45x10^ 

Using  a  three-parameter  (Eq,  Ei  ,  and  Xj )  viscoelastic  model,  the  relaxation  function  takes  the 
form 


T(t)=Eo+Eie'''^' 

(9) 

The  equivalent  loss  factor,  T),  for  a  viscoelastic  material  with  the  relaxation  modulus  shown 
above  is  as  follows: 

E" 

E' 

(10) 

where  E'  and  E"  are  the  storage  and  loss  moduli,  respectively,  and  can 

be  found  as  follows: 

E'  =  E,+(c»T|y  ,  ' 

(1+(®T,)  j 

(11) 

(12) 

with  0)  equal  to  the  loading  frequency  in  rad/s.  Earlier  work  [12]  indicates  that  the  dissipation 
due  to  the  deviatoric  and  dilatational  strains  are  of  similar  magnitudes.  As  a  result,  the 
relaxation  function  is  not  broken  into  separate  deviatoric  and  dilatational  components  or, 
equivalently,  it  is  assumed  that  the  Poisson’s  ratio  remains  constant. 

Figure  30  shows  elastic  modulus  and  loss  factor  versus  frequency  for  a  viscoelastic  material 
with  Eo  equal  to  24.0x10^  psi,  Ei  equal  to  12.0x10®  psi,  and  Ti  equal  to  2.0x10®.  These 
properties  correspond  to  a  maximum  material  loss  factor  of  approximately  0.200  and  an  elastic 
modulus  of  approximately  30.0x10®  psi  at  the  frequency  of  the  maximum  loss  factor.  (Note  that 
these  properties  are  not  representative  of  the  typical  intrinsic  damping  of  steel,  which  may  be 
much  less  than  0.200,  but  have  been  used  for  illustrative  purposes.)  The  material  has  an 
instantaneous  modulus  of  36.0x10®  psi  and  a  long-term  modulus  of  24.0x10®  psi.  Based  on  the 
results  shown  inTable  3,  it  is  likely  that  the  majority  of  the  impacts  during  a  particle  damper 
simulation  will  occur  at  frequencies  between  approximately  1.45x10"^  and  1.45x10  hertz.  This 
range  is  shown  in  Figure  30. 
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Figure  30.  Modulus  and  Loss  Factor  versus  Frequency  for  Three-Parameter  Model 

A  five-parameter  (Eo,  Ei,  Ti,  E2,  and  T2),  or  higher  parameter,  viscoelastic  model  could  also  be 
used  to  simulate  the  viscoelastic  behavior  at  the  expense  of  additional  equations  which  must  be 
solved  at  every  contact  increment.  Figure  31  shows  elastic  modulus  and  loss  factor  versus 
frequency  for  a  viscoelastic  material  with  Eo  equal  to  20.0x10^  psi,  Ei  equal  to  12.0x10  psi,  Xi 
equal  to  1.0x10'®,  E2  equal  to  8.0x10®  psi,  and  t2  equal  to  1.0x10'®.  These  properties  also 
correspond  to  a  maximum  material  loss  factor  of  approximately  0.200  and  an  elastic  modulus 
of  approximately  30.0x10®  psi  at  the  frequency  of  the  maximum  loss  factor.  A  three-parameter 
representation  has  been  incorporated  into  the  current  particle  damper  model;  however,  the 
model  could  easily  be  extended  to  incorporate  higher  parameter  representations  at  the  expense 
of  additional  computational  resources. 
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Figure  31.  Modulus  and  Loss  Factor  versus  Frequency  for  Five-Parameter  Model 

Using  the  existing  finite  element  model,  ABAQUS  analyses  have  been  performed  to  assess  the 
effect  of  using  viscoelastic  properties  for  the  spheres.  The  steel  sphere  has  been  given  the  three- 
parameter  viscoelastic  properties  as  shown  in  Figure  30.  Figure  32  shows  predicted  load  and 
displacement  as  a  function  of  time  when  elastic  properties  are  used  and  when  viscoelastic 
properties  are  used.  Note  that  when  viscoelastic  properties  are  used,  a  lesser  maximum  total 
force  is  predicted,  with  the  peak  force  occurring  prior  to  the  peak  displacement.  After  the  peak 
displacement,  the  displacements  tend  to  remain  larger  longer  and  require  a  longer  time  to  return 
to  zero. 
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Figure  32.  Predicted  Elastic  and  Viscoelastic  Behavior 

It  is  desirable  to  express  the  force-displacement  relations  for  the  dissipative  portion  of  the 
normal  force  in  similar  form  as  the  relations  for  the  elastic  portion,  shown  previously  in 
Equation  5.  The  force-displacement  relations  used  by  Saluena,  et  al.,  have  been  based  on  the 
following  expression  derived  by  several  researchers  [15-16]  by  assuming  viscoelastic 
properties  for  the  spheres 

F5.)  =  A^R'V”d  (13) 

In  Equation  13,  A  is  a  material  constant  which  characterizes  the  dissipative  character  of  the 
material.  In  general,  however,  viscoelasticity  is  a  function  of  time  (or  frequency).  A  is  given  as 
a  function  of  various  parameters  describing  the  viscoelastic  behavior  of  the  particles  (more 
precisely,  the  coefficients  of  viscosity  relating  to  the  deviatoric  and  dilatational  strains).  These 
parameters  may  be  difficult  to  obtain  and,  as  a  result,  previous  efforts  have  used  A  as  a  fit 
parameter  [15].  To  predict  damping,  it  will  be  necessary  to  define  the  viscoelastic  properties  of 
the  sphere  material  prior  to  performing  the  mathematical  simulation. 

Ideally,  the  dissipative  portion  of  the  normal  force  would  be  expressed  in  terms  of  a 
viscoelastic  formulation  of  Hertz’s  theory.  Recall  that  Hertz’s  theory  is  based  on  the  basic 
equations  of  elasticity  and  the  geometries  of  the  contact  bodies.  When  linear  isotropic 
viscoelasticity  is  assumed,  only  the  constituitive  relations  of  the  basic  elasticity  equations  are 
changed.  The  constituitive  relations  can  be  expressed  in  general  form  as  follows: 

Pa  =  Qe  (14) 
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For  viscoelasticity,  P  and  Q  are  linear  operators  in  the  time  variable  which  express  the 
viscoelastic  behavior  of  the  material.  These  operators  may  include  differential  operators, 
integral  operators,  or  any  other  equivalent  means  of  expressing  viscoelastic  behavior.  For  the 
three-parameter  viscoelastic  model  used  for  this  effort,  the  following  applies: 

P  =  1  (15) 

Q  =  T(t)=Eo-hE,e'‘'^‘  (16) 

and  the  constituitive  relation  (relating  the  normal  stress  and  normal  strain)  can  be  written  as 
follows: 

<7  =  'r(j>  =  ^„+E,e-''''>  (17) 

The  assumptions  inherent  in  the  Hertz  solution  also  must  not  be  violated.  In  particular,  these 
assumptions  must  ensure  that:  (1)  the  boundaries  of  the  bodies  do  not  significantly  affect  the 
stresses  in  the  contact  zone;  (2)  the  portions  of  the  bodies  outside  of  the  contact  zone  can  be 
roughly  approximated  as  a  viscoelastic  half-space;  and,  (3)  the  strains  remain  sufficiently  small 
that  the  linear  theory  of  viscoelasticity  applies.  To  ensure  the  above  conditions,  the  dimensions 
of  the  contact  area  must  be  small  compared  with  the  dimensions  of  each  body  and  with  the 
relative  radii  of  curvature  of  the  surfaces. 

Pao  [17]  was  one  of  the  first  to  attempt  to  extend  the  Hertz  solution  to  include  linear 
viscoelasticity.  Pao  took  the  Laplace  transform  of  the  viscoelastic  equations  (i.e.,  the  basic 
elasticity  equations  with  a  viscoelastic  constituitive  equation),  to  remove  the  time  variable  and 
yield  an  elastic  problem  in  terms  of  transformed  variables.  This  solution  method  is  valid  if  the 
regions  over  which  different  types  of  boundary  conditions  are  prescribed  do  not  vary  during  the 
time  under  consideration.  However,  for  the  special  case  of  the  impact  of  viscoelastic  spheres, 
these  regions  vary  and  this  approach  fails  since  for  some  points  on  the  surface,  neither  the 
traction  nor  the  displacement  is  known  throughout  the  history  of  the  problem,  so  the  transform 
of  neither  can  be  determined  [18]. 

Radok  [19]  suggested  an  alternative  approach  of  taking  a  one-parameter  family  of  solutions  of 
the  elastic  problem  in  the  time  domain,  with  the  same  boundary  conditions  as  the  viscoelastic 
problem,  and  replacing  the  elastic  constants  with  the  appropriate  viscoelastic  operators  in  the 
expressions  for  the  stress  components.  This  procedure  can  permit  the  manipulation  of  the 
elastic  constants  in  the  determination  of  the  elastic  solution  to  involve  procedures  which  are  not 
valid  for  operators,  or  which  lead  to  the  introduction  of  new  functions  in  the  operator  case.  As  a 
result,  the  significance  of  any  solution  must  be  checked,  often  by  changing  the  initial 
formulation  so  that  a  direct  transform  attack  is  available.  Using  this  alternative  approach,  Lee 
and  Radok  [18]  have  shown  that,  as  long  as  the  contact  area  is  increasing,  a  simple  relation  for 
the  normal  force  can  be  derived  by  replacing  the  elastic  modulus  in  Hertz’s  relation  (Equation 
5)  with  the  relaxation  function  for  the  sphere  material.  Substituting  the  relaxation  function  into 
Hertz’s  relation  and  recognizing  that  the  contact  radius  also  is  a  function  of  time,  the  total 
normal  force  (i.e.,  the  normal  force  due  to  the  combined  elastic  and  dissipative  components)  at 
any  time  can  be  expressed  as  follows: 
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Substituting  the  relation  in  Equation  7  for  the  contact  radius,  the  derivative  in  Equation  18  can 
be  evaluated.  Equation  9  also  can  be  substituted  for  the  relaxation  function  yielding,  as  follows; 

'p  1/2  t  n 

FMt)=T^jiE,e-(‘-‘>^a(t7'^d(t>'  (19) 

(l-t)  j;'i=o 

For  greater  utility,  it  would  be  beneficial  to  express  Equation  19  in  incremental  form  such  that 
the  total  load  at  any  time  can  be  calculated  based  only  on  information  from  the  previous  time. 
Otherwise,  it  will  be  necessary  to  retain  the  entire  loading  history.  The  contact  force  at  a  time 
equal  to  the  original  time,  t,  plus  an  additional  increment  in  time.  At,  can  be  expressed  as 
follows: 


F^(t  +  At)= 


n 

J  jE.e-(‘""'-‘>‘a(t'f  d(t')dt' 


(20) 


The  integral  in  Equation  20  can  be  broken  into  two  parts,  one  from  time  0  to  t,  and  the  second 
from  time  t  to  t+At,  as  follows: 


F^(t  +  At)  = 


a(t'rd(t>'  + 


0  i=0 


(21) 


By  considering  the  individual  contributions  of  each  of  the  viscoelastic  terms  to  the  total  contact 
force,  the  total  load  at  any  time  can  be  expressed  as  a  summation  of  terms,  as  follows: 

F^(t)=  jFi'^(t)  (22) 

i=0 


Further,  the  first  term  in  brackets  on  the  right  side  of  Equation  21  is  the  contact  force  at  time  t 
(as  shown  in  Equation  19)  multiplied  by  a  decaying  exponential.  As  a  result,  each  of  the 
individual  terms  contributing  to  the  contact  force  can  be  expressed  as  follows: 


Fi''(t  +  At)=Fi''(t>-^‘'’‘ 


p  1/2  t+At 

.  f a(t')''' d(t')dt' 
(1-u  j  J 


+ 


(23) 


Using  the  midpoint  rule  to  evaluate  the  integral  on  the  right-hand  side  of  Equation  23  yields  the 
following: 


Fi^(t  +  At)=Fi''(t> 


-At/Xj  ^ 


1/2 


For  the  special  case  of  i=0,  the  following  applies: 


Eie-^'^^'a(t  +  At/2)‘''d(t  +  At/2)^t  (24) 


Fj^(t  +  At)=F(r(t> 


R 


1/2 


+  At/2)\t 


(25) 


For  other  values  of  i: 
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Fi""  (t  +  At)=  +  At/2)'''d(t  +  At/2)At  (26) 


Equations  25  and  26  can  be  used  to  find  the  contact  force  provided  the  contact  area  is 
increasing.  However,  when  the  loading  history  causes  the  contact  area  to  decrease,  this 
technique  breaks  down.  The  breakdown  results  due  to  unrealistic  negative  contact  pressures 
which  are  predicted  within  the  contact  area.  Ting  [20-21]  studied  the  case  when  the  contact  area 
is  decreasing  and  found  that  the  impact  could  be  described  by  the  following  pair  of  coupled 
integro-differential  equations  for  the  approach  and  the  total  contact  load: 


Ra(t)=a(t)"-Jo(t-t')^ 


jT(t'-r)A^(tO^> 


dt' 


(27) 


In  Equations  27  and  28,  t  is  the  actual  time  at  which  the  equations  are  being  evaluated,  t*  is  the 
time  of  maximum  penetration,  and  ti  is  time  during  the  loading  (i.e.,  ti<t )  at  which  the  contact 
size,  a(ti),  was  equal  to  the  current  contact  size,  a(t).  As  previously  discussed,  0(t)  is  the  creep 
compliance  and  ^(1)  is  the  relaxation  function.  Because  the  contact  size  at  time  ti  is  required. 
Equations  33  and  34  must  be  evaluated  numerically  step  by  step  [22].  Since  the  entire  loading 
history  must  be  stored  and  the  number  of  required  function  evaluations  is  constantly  increasing 
during  the  unloading,  this  technique  requires  considerable  memory  and  computational 
resources.  Figure  33  shows  load  and  displacements  curves  as  predicted  using  Equations  27  and 
28  along  with  finite  element  analysis  results  for  comparison. 


40 


Figure  33.  Calculated  Viscoelastic  Behavior  Using  Ting’s  Method 

As  seen  in  Figure  33,  there  is  a  fairly  significant  difference  between  the  predicted  results  using 
Ting’s  method  and  those  from  the  ABAQUS  finite  element  analysis,  particularly  in  the  loads 
late  in  the  impact  event.  Using  the  procedure  outlined  by  Calvitt  [22]  to  solve  Equations  27  and 
28,  the  force  and  displacement  are  predicted  to  return  to  zero  at  the  same  time.  Due  to  the 
relaxation  of  the  viscoelastic  material,  one  would  expect  the  force  to  return  to  zero  prior  to 
reaching  zero  displacement  (i.e.,  the  particle  and  cavity  would  still  be  in  contact  although  no 
force  would  be  transmitted  through  the  contact). 

An  alternative  approach  which  greatly  reduces  the  required  resources  is  to  calculate  the  total 
contact  force  using  Equations  25  and  26,  with  the  assumption  that  the  total  contact  force  goes 
to  zero  when  negative  forces  are  predicted.  As  shown  in  the  load  and  displacements  curves  in 
Figure  34,  the  viscoelastic  behavior  calculated  using  this  method  compares  well  with  that 
predicted  using  finite  element  analysis  and  the  predicted  force  and  displacement  do  not  return 
to  zero  at  the  same  time.  This  method  has  been  incorporated  in  the  particle  damper  model.  Note 
that  these  relations  capture  both  the  elastic  and  dissipative  portions  of  the  normal  force  and 
revert  back  to  an  incremental  form  of  Hertz’s  relation  when  elastic  properties  are  used. 
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Figure  34.  Calculated  Viscoelastic  Behavior  Using  Modified  Radok’s  Method 
7.2.4  Shear  Force 

In  addition  to  normal  forces,  shear  forces  also  can  be  generated  between  individual  particles  in 
contact  or  particles  in  contact  with  the  cavity  walls.  In  the  absence  of  friction,  and  for  the 
impact  of  two  bodies  with  identical  material  properties,  the  normal  and  shear  forces  are 
independent  from  each  other.  When  friction  is  introduced,  or  the  impacting  bodies  are  no 
longer  of  the  identical  material,  the  normal  and  shear  forces  interact.  However,  this  interaction 
is  generally  small  and  it  typically  is  assumed  that  the  normal  and  shear  forces  are  independent 
of  each  other  [23].  This  assumption  is  made  in  the  particle  damper  model. 

Shear  forces  are  created  when  a  particle  slides  along  another  particle  or  along  the  cavity  walls, 
or  when  oblique  impacts  occur  between  individual  particles  or  between  particles  and  the  cavity 
walls.  Both  slipping  and  nonslipping  contact  can  occur.  The  enduring  contacts  seen  as  a  particle 
slides  along  another  particle  or  along  the  cavity  walls  typically  last  significantly  longer  in  time 
than  the  nearly  instantaneous  oblique  impacts  between  individual  particles  or  between  particles 
and  the  cavity  walls.  As  a  result,  the  enduring  contacts  can  have  a  much  more  significant  effect 
on  the  overall  damping  than  the  oblique  impacts,  and  it  is  important  that  the  model  accurately 
capture  these  effects. 

One  of  the  simplest  models  which  incorporates  slipping  at  the  point  of  contact  is  Amonton’s 
law  of  sliding  friction,  or  Coulomb  friction  [23]: 

F'=-sgn(v:„V|F''|  (29) 
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where  |a  represents  the  particle-particle  friction  coefficient  and  is  the  relative  tangential 
velocity.  This  simple  model  generally  performs  well;  however,  problems  can  arise  with  this 
model  when  the  slip  rate  is  small,  as  the  shear  force  can  rapidly  change  directions  and  create 
oscillations.  A  second  problem  with  this  simple  model  is  the  inability  to  model  static  friction 
[24]. 


To  avoid  force  oscillations  possible  in  the  Coulomb  friction  model,  previous  researchers, 
including  Saluena,  et  al.  [10],  have  included  a  “viscous  friction”  in  the  formulation.  Haff  and 
Werner  [25]  proposed  the  following  force-displacement  relation: 


=  -sgn 


(v;„)ninf 


Y.meff 


'rel 


(30) 


with  the  effective  mass,  m^ff ,  found  as  follows: 


K  +m/) 


(31) 


where  nij  and  nij  are  the  masses  of  the  contacting  bodies.  The  niin{  }  term  in  Equation  26 

ensures  that  the  maximum  no-slip  friction  force  cannot  exceed  the  Coulomb  friction  force  as 
given  in  Equation  25.  A  shear  damping  coefficient,  Vt .  is  used  to  introduce  a  “viscous  friction” 
force  proportional  to  the  relative  velocity  and  the  effective  mass.  This  ensures  that  the 
tangential  force  equation  is  continuous  and  prevents  the  possibility  of  large  force  oscillations  if 
only  the  Coulomb  friction  force  were  used.  However,  this  model  also  cannot  simulate  static 
friction  [24]  and  requires  the  definition  of  the  shear  damping  coefficient  which  artificially  adds 
damping  to  the  system. 

Herrman  [26]  proposed  the  following  relation  incorporating  a  tangential  virtual  spring: 

=-sgn(5')nin^k‘5‘, pip'll]  (32) 


where  k*  is  the  stiffness  of  the  tangential  virtual  spring  and  5‘  is  the  total  shear  displacement 
during  the  contact  and  is  found  as  follows: 

=  J  vj,,dt  (33) 

In  this  method,  when  two  particles  initially  contact,  the  virtual  spring  is  activated  in  the 
tangential  direction.  This  spring  creates  a  restoring  force  until  the  tangential  force  exceeds  the 
force  given  by  the  Coulomb  criteria,  after  which  the  spring  is  removed  and  sliding  occurs.  A 
similar  procedure,  typically  referred  to  as  elastic  slip,  is  used  in  a  number  of  commercial  finite 
element  codes  (such  as  ABAQUS  [27])  when  solving  contact  problems  with  friction  present. 
This  model  captures  static  friction,  but  oscillations  again  are  possible  if  the  relative  tangential 
displacement  is  in  the  opposite  direction  as  the  relative  tangential  velocity  [24].  The  selection 
of  an  appropriate  value  for  the  stiffness  of  the  tangential  virtual  spring  also  is  a  concern. 

All  of  these  models  revert  back  to  the  Coulomb  friction  model  when  significant  sliding  occurs. 
As  a  result,  the  Coulomb  friction  model  has  been  incorporated  into  the  particle  damper  model 
to  account  for  the  shear  forces  created  during  enduring  contacts.  A  single  coefficient  of  friction 
is  used  for  particle-particle  and  particle-cavity  contact.  The  kinetic  coefficient  of  friction  is 
given,  and  it  is  assumed  that  relative  particle  motion  occurs  or  that  the  static  coefficient  of 
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friction  is  zero.  This  assumption  is  reasonable  since  relative  particle  motion  must  occur  for  the 
particle  damper  to  function. 

Under  gravity  loads  or  other  body  loads  (e.g.,  centrifugal  loads),  it  is  possible  that  particle 
motion  may  be  inhibited.  As  a  result,  it  is  recommended  that  some  simple  first-order 
calculations  be  performed  to  determine  if  particle  motion  is  likely.  Vibratory  forces  on  a 
particle  can  be  estimated  based  on  the  mass  of  the  particle  and  the  vibratory  accelerations. 
Friction  forces  can  be  estimated  using  the  Coulomb  friction  model  with  the  coefficient  of 
friction  defined  and  the  normal  forces  estimated  based  on  the  mass  of  the  particle  and  the 
gravitational  (or  centrifugal)  acceleration.  Comparison  of  these  forces  (i.e.,  ensuring  that  the 
vibratory  forces  are  sufficient  to  overcome  the  frictional  forces)  should  provide  a  good 
indication  as  to  the  likelihood  of  particle  motion. 

To  investigate  the  effects  of  oblique  impacts,  finite  element  analyses  of  an  elastic  sphere 
impacting  a  rigid  plate  have  been  performed.  Due  to  the  unsymmetric  loads  generated  during 
the  impact,  a  three-dimensional  half-symmetry  model,  shown  in  Figure  35,  has  been  used  for 
the  analyses.  The  model  contained  37,138  elements  (C3D8  elements  for  the  sphere  and  R3D4 
elements  for  the  rigid  plate)  to  accurately  capture  the  contact  behavior.  Preprocessing  and 
postprocessing  of  the  finite  element  model  and  results  was  performed  using  MSC/PATRAN 
Version  7.6  from  the  MacNeal  Schwendler  Corporation  [28]  with  the  analyses  performed  using 
Version  5.7  of  the  ABAQUS/Standard  finite  element  code  from  Hibbitt,  Karlsson,  and 
Sorensen,  Incorporated  [27]. 


Figure  35.  Three-Dimensional  Model  of  Oblique  Impact  of  Spherical  Particle  and  Rigid 

Wall 
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Initial  analyses  investigated  an  elastic  sphere  of  radius  0.0625  inch  impacting  the  rigid  wall  at  a 
velocity  of  62.832  in/s  at  an  angle  of  45°  relative  to  the  rigid  wall.  The  sphere  was  given  steel 
properties  with  an  elastic  modulus  of  30x10^  Ibf/in^,  a  density  of  0.283  Ibf/in^,  and  a  Poisson  s 
ratio  of  0.33.  A  coefficient  of  friction  of  0.60  was  prescribed.  Predicted  normal  and  shear  loads 
resulting  due  to  the  oblique  impact  are  shown  in  Figure  36. 


Time  (sec) 


Figure  36.  Predicted  Normal  and  Shear  Loads  Resulting  from  Oblique  Impact 

Figure  37,  Figure  38,  and  Figure  39  show  predicted  normal  and  shear  loads  versus  time  from 
ABAQUS  and  as  predicted  using  the  Coulomb  friction  model  for  impacts  at  angles  of  incidence 
of  60°,  45°,  and  15°,  respectively  (where  the  angle  of  incidence  is  the  angle  between  the 
particle  velocity  vector  and  the  direction  perpendicular  to  the  rigid  wall).  The  normal  forces 
were  calculated  based  on  Hertz’s  theory  and,  as  previously  discussed,  compare  well  with  the 
results  from  ABAQUS.  As  shown  in  the  figures,  the  Coulomb  friction  model  generally  predicts 
reasonably  accurate  shear  forces  for  impacts  at  higher  angles  of  incidence  (more  glancing). 
However,  at  lower  angles  of  incidence  (more  normal)  the  Coulomb  friction  model  can 
significantly  overpredict  the  shear  forces.  In  addition,  the  shape  of  the  shear  force  curve  from 
ABAQUS  is  different  than  that  predicted  using  the  Coulomb  friction  model.  In  the  ABAQUS 
solution,  the  shear  force  peaks  prior  to  the  peak  of  the  normal  force.  With  the  Coulomb  friction 
model,  the  shear  and  normal  forces  peak  at  the  same  time.  The  shear  force  also  undergoes  a 
reversal  in  the  ABAQUS  solution,  whereas  the  normal  force  completes  a  half  cycle  only  [23]. 
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re  37.  Predicted  Normal  and  Shear  Loads  Using  Coulomb  Friction  for  60°  Impact 


Figure  38.  Predicted  Normal  and  Shear  Loads  Using  Coulomb  Friction  for  45°  Impact 
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Figure  39.  Predicted  Normal  and  Shear  Loads  Using  Coulomb  Friction  for  15°  Impact 

Investigating  the  shear  forces  further,  one  finds  that  the  shear  force-displacement  relations  are 
extremely  complex.  For  example,  the  force-displacement  relations  can  independently  display 
nonlinear  response,  path  dependence,  and  dissipation  due  to  slip,  even  when  the  spheres  are 
linearly  elastic  and  the  relative  displacements  are  small.  Mindlin  and  Deresiewicz  [29] 
established  a  detailed  set  of  equations  depending  on  whether  the  normal  and  shear  forces  are 
increasing,  decreasing,  or  remaining  constant.  Their  theory  is  limited  to  so-called  simple 
loading  histories  and  includes  11  loading  cases.  Simple  loading  histories  begin  with  an 
equilibrium  position  for  which  the  distribution  of  traction  is  equivalent  to  a  state  of  constant 
normal  force  and  varying  shear  force,  and  advance  to  the  desired  state  through  a  sequence  of 
equilibrium  positions.  The  final  distribution  of  traction  is  equivalent  to  another  state  of  constant 
normal  force  and  varying  shear  force.  Due  to  the  complexity  of  these  relations,  efficient 
implementation  within  the  particle  dynamics  method  is  not  possible.  As  a  result,  various 
simplified  contact  force-displacement  relations  have  been  proposed. 

Walton  and  Braun  [30]  proposed  an  incremental  equation,  based  on  a  tangential  stiffness 
coefficient  and  the  incremental  tangential  displacement,  to  approximate  the  Mindlin  and 
Deresiewicz  contact  mechanics  theory.  For  cases  where  the  normal  force  remains  constant,  this 
method  generates  results  in  general  agreement  with  those  from  Mindlin  and  Deresiewicz. 
However,  for  particle  damping  simulations,  and  the  majority  of  particle  dynamics  simulations, 
both  the  normal  and  shear  forces  are  changing  during  the  impact  event. 

Vu-Quoc  and  Zhang  [31-33]  have  proposed  an  improved  version  of  the  model  by  Walton  and 
Braun  which  incorporates  selected  cases  of  varying  normal  and  shear  forces  from  Mindlin  and 
Deresiewicz,  and  a  correction  for  the  rolling  effect.  Although  these  relations  are  simpler  than 
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those  of  Mindlin  and  Deresiewicz,  considerable  computational  resources  still  are  required.  As 
an  example,  Vu-Quoc  and  Zhang  [33]  performed  simulations  of  400  soybeans  falling  through  a 
chute,  which  required  approximately  19  CPU  hours  for  1  second  of  simulated  time.  As 
discussed  in  the  following  Experimental  Validation  and  Sample  Application  sections,  particle 
damper  simulations  may  require  solutions  for  tens  of  seconds.  The  computational  resource 
requirements  make  implementation  of  the  Vu-Quoc  and  Zhang  relations  impractical. 

A  simpler,  more  efficient  method  which  can  reasonably  account  for  the  tangential  force- 
displacement  relations  is  required.  Maw,  Barber,  and  Fawcett  [34-35]  proposed  a  solution  for 
the  oblique  impact  of  an  elastic  sphere  where  the  contact  area  is  broken  into  a  series  of 
concentric  circles.  In  each  of  these  separate  regions,  sticking  or  slipping  can  occur.  At  the 
beginning  of  each  step  in  the  solution,  each  of  the  regions  is  assumed  to  either  stick  or  slip.  The 
appropriate  equations  are  solved  and  the  solution  is  tested  to  see  if  the  initial  assumptions  were 
correct.  In  stick  regions,  the  tangential  force  must  not  exceed  the  limiting  force  given  by  the 
Coulomb  friction  model.  In  slip  regions,  the  incremental  displacement  must  be  in  the  correct 
sense  for  the  assumed  frictional  traction.  If  either  of  these  tests  fails  in  any  region,  the  initial 
assumptions  are  revised  and  a  new  solution  is  obtained.  At  each  contact  increment,  a  solution 
must  be  found  for  each  of  the  different  regions  (which  may  be  on  the  order  of  20);  each 
solution  must  be  checked  to  verify  that  the  initial  assumptions  for  that  particular  region  are 
correct;  and,  as  necessary,  new  solutions  must  be  found  for  those  regions  in  which  the  initial 
assumptions  were  violated.  As  a  result,  this  method  is  relatively  inefficient. 

Maw,  Barber,  and  Fawcett  extended  their  work  by  calculating  the  shear  forces  for  various 
oblique  impact  conditions  and  found  that  the  variation  in  the  shear  force  during  the  impact 
could  be  found  based  on  two  nondimensional  parameters.  These  efforts  are  outlined  in 
Johnson’s  Contact  Mechanics  [23]  and  have  the  advantage  of  using  the  non-dimensional 
parameters  to  turn  the  computationally  intensive  procedure  into  one  where  results  can  be  found 
to  fall  on  a  specific  curve.  The  peak  shear  forces  predicted  using  the  Maw,  Barber,  and  Fawcett 
model  compare  relatively  well  with  the  peak  shear  forces  from  ABAQUS  and  the  model 
captures  the  basic  form  of  the  curve  including  the  reversal.  The  peak  shear  force  predicted 
using  the  Maw,  Barber,  and  Fawcett  model  occurs  earlier  in  the  impact  event  than  that 
predicted  by  ABAQUS  and,  particularly  for  intermediate  angles  of  incidence,  a  much  larger 
negative  shear  force  is  predicted  during  the  reversal. 

All  of  the  shear  force-displacement  models  discussed  above  suffer  from  the  same  general  flaw 
that,  in  general,  it  is  not  possible  to  know  a  priori  whether  a  contact  will  be  enduring  or  not. 
The  relations  discussed  above  are  based  on  a  single  discrete  impact,  and  may  not  accurately 
predict  the  shear  forces  for  enduring  contacts  or  for  particles  experiencing  multiple  contacts  at 
once.  In  the  particle  damper  model,  it  is  more  critical  to  accurately  represent  the  enduring 
contacts  since  these  contacts  generally  last  for  a  longer  time  and  can  have  a  more  significant 
effect  on  the  damping.  As  a  result,  shear  force-displacement  relations  for  oblique  impacts  have 
been  implemented  in  the  particle  damping  model  using  the  same  Coulomb  friction  model  used 
for  enduring  contacts.  Future  enhancements  to  the  particle  damper  model  may  include  a  more 
accurate  representation  of  the  shear  forces  generated  during  oblique  impacts. 

7.3  Implementation 

The  background  of  the  particle  damper  simulation  code  is  based  on  X3D,  an  explicit  finite 
element  code  typically  used  for  impact  analyses  [36].  The  code  contains  various  contact 
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algorithms  and  bookkeeping  routines  and  provides  an  appropriate  framework  for  simulating 
particle  damping  through  the  use  of  the  particle  dynamics  method.  Particle-particle  and 
particle-cavity  contacts  are  resolved  using  the  force-displacement  relations  discussed  in  the 
preceding  paragraphs.  Separate  contact  detection/resolution  routines  are  used  for  particle- 
particle  and  particle-cavity  contacts. 

Particle-cavity  contact  detection  is  performed  using  a  modified  version  of  the  master-slave 
algorithm  from  X3D.  As  shown  in  Figure  40,  three-node  triangular  master  surfaces  are  used  to 
define  the  cavity  with  the  particles  existing  as  slave  nodes.  At  each  time  step,  geometry  checks 
are  performed  to  resolve  any  contact  between  the  slave  nodes  (particles)  and  the  master 
surfaces  (cavity).  Since  the  cavity  walls  are  assumed  to  be  flat,  the  nodes  used  for  the  master 
surface  lie  in  a  plane.  An  orthogonal  local  coordinate  system  is  constructed  with  two  directions, 
Xi  and  X2,  as  shown  in  the  figure,  in  the  plane  of  the  master  surface  and  the  third  direction,  X3, 
perpendicular  to  the  contact  surface.  The  distance  in  the  X3  (normal)  direction  between  the 
particle  center  and  the  master  surface  is  calculated.  Contact  is  detected  if  this  distance  is  less 
than  the  particle  radius.  The  particle  translational  velocity,  Vp,  is  broken  into  components  in  the 
local  coordinate  system.  Similarly,  the  master  node  velocities  (v^i,  Vm2,  and  Vms)  are  broken 
into  translation  components  in  the  local  coordinate  system.  Any  particle  rotation,  (Op,  is  broken 
into  rotations  about  the  in-plane  (xi  and  X2)  axes.  Spin  (i.e.,  rotation  about  the  X3  axis)  is 
neglected.  The  relative  normal  velocity  is  calculated  based  on  the  X3  component  of  the  particle 
velocity  and  X3  contributions  from  each  of  the  master  nodes.  Relative  tangential  velocities  are 
calculated  based  on  the  in-plane  velocity  components  of  the  particle  and  the  in-plane  velocity 
contributions  from  each  of  the  master  nodes,  as  well  as  the  particle  rotations  about  the  in-plane 
axes.  A  slight  error  is  introduced  as  the  velocity  contributions  from  the  particle  rotation  about 
the  in-plane  axes  are  based  on  undeformed  particle  radii  rather  than  the  instantaneous  radii. 


Figure  40.  Particle-Cavity  Contact  Detection  and  Resolution 

Particle-particle  contact  detection  is  performed  by  tracking  the  nodes  at  the  center  of  the 
particles  as  shown  in  Figure  41.  The  distance  between  each  set  of  particle  centers  is  calculated. 
Contact  is  detected  if  this  distance  is  less  than  twice  the  particle  radius.  The  relative  velocity  at 
the  point  of  contact  is  calculated  in  global  coordinates  based  on  the  translations  (vpi  and  Vp2) 
and  rotations  (tOpi  and  (Op2)  of  the  contacting  particles.  As  with  the  particle-cavity  contacts,  a 
slight  error  is  introduced  as  the  velocity  contributions  from  the  particle  rotations  are  based  on 
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undeformed  particle  radii  rather  than  the  instantaneous  radii.  A  unit  normal  vector  in  the 
direction  connecting  the  particle  centers  is  calculated  in  global  coordinates.  The  relative  normal 
velocity  can  be  calculated  by  taking  the  dot  product  of  the  relative  velocity  vector  at  the  point 
of  contact  with  the  unit  normal  vector.  The  relative  tangential  velocity  then  can  be  calculated 
by  subtracting  the  relative  normal  velocity  from  the  overall  relative  velocity. 


Figure  41.  Particle-Particle  Contact  Detection  and  Resolution 

The  force-displacement  relations  used  to  resolve  both  particle-cavity  and  particle-particle 
contacts  are  formulated  in  incremental  form.  As  a  result,  whenever  contact  is  detected,  it  is 
necessary  to  determine  whether  it  is  a  new  or  existing  contact  condition.  Checks  are  made  to 
determine  if  the  contact  condition  existed  at  the  previous  time  step  and  appropriate  information 
from  the  previous  time  step,  which  is  required  to  calculate  the  forces  at  the  current  time  step,  is 
stored.  A  constant  time  step  has  been  used  throughout  the  simulations. 

7.4  Summary  of  Model  Assumptions  and  Limitations 

The  following  is  a  summary  of  the  primary  assumptions  and  limitations  of  the  particle  damper 
model: 

•  All  particles  assumed  to  be  identical;  particles  are  assumed  to  be  spherical  in  shape  with 
a  given  radius  and  material  properties. 

•  Particles  can  be  given  elastic  or  viscoelastic  material  properties;  viscoelastic  material 
properties  are  given  as  a  three-parameter  Maxwell  model. 

•  Cavity  walls  are  assumed  to  be  flat  and  rigid. 

•  A  single  coefficient  of  friction  is  used  for  particle-particle  and  particle-cavity  contact. 

•  The  kinetic  coefficient  of  friction  is  given;  i.e.,  it  is  assumed  that  particle  motion  occurs 
or  that  the  static  coefficient  of  friction  is  zero. 

•  Normal  force-displacement  relations  are  based  on  a  modified  form  of  Lee  &  Radok’s 
relation  (where  the  elastic  modulus  is  replaced  with  the  relaxation  modulus)  with  the 
force  set  to  zero  when  a  negative  force  is  predicted. 

•  For  elastic  properties,  the  normal  force-displacement  relation  reverts  to  an  incremental 
form  of  Hertz’s  relation. 


•  For  viscoelastic  properties,  an  incremental  form  of  the  force-displacement  relation  is 
used  which  incorporates  relaxation  behavior  of  the  viscoelastic  material. 

•  Shear  force-displacement  relations  are  based  on  Amonton’s  law  of  sliding  friction,  or 
Coulomb  friction. 

•  The  shear  force  magnitude  is  based  solely  on  the  magnitude  of  the  normal  force  and  the 
coefficient  of  friction. 

•  The  direction  of  the  shear  force  opposes  the  relative  tangential  velocity  between  the 
contacting  surfaces;  relative  tangential  velocities  can  result  from  oblique  impacts  or  due 
to  rotation  of  the  particles. 

7.5  Particle  Damper  Design  Methodology 

As  discussed  in  the  preceding  paragraphs,  an  analytical  particle  damper  model  has  been 
developed  which  captures  the  complex  interactions  of  the  loss  mechanisms  in  a  particle 
damper.  This  simulation  technique  has  been  incorporated  into  a  comprehensive  particle  damper 
design  methodology  which  allows  particle  damping  to  be  implemented  without  extensive  trial- 
and-error  testing.  Critical  steps  in  the  design  process  include  the  following: 

1.  Determine  characteristics  of  undamped  system. 

2.  Determine  appropriateness  of  particle  dampers. 

3.  Select  particle  damper  configuration. 

4.  Determine  characteristics  of  undamped  system  with  adjusted  mass. 

5.  Evaluate  damper  effectiveness  using  particle  damper  simulation  technique. 

6.  As  necessary,  repeat  steps  3-6. 

7.  Verify  chosen  design  experimentally. 

The  initial  step  in  the  particle  damper  design  is  to  establish  the  charaeteristics  of  the  system  to 
be  damped.  Key  characteristics  include  the  resonant  frequencies  of  the  structure  and  the 
vibratory  displacements  of  the  system.  These  characteristics  will  typically  be  determined 
through  a  combination  of  experimental  testing  and  finite  element  analysis.  Additional  concerns 
include  the  environmental  conditions  under  which  the  damper  must  operate  (e.g.,  at  cryogenic 
or  elevated  temperatures).  Geometric  constraints,  which  may  limit  the  location  and/or  size  of 
the  particle  damper  cavity,  also  must  be  considered. 

The  second  step  in  the  design  methodology  is  to  determine  if  particle  dampers  are  appropriate. 
Particle  damping  is  typically  most  effective  in  damping  well-separated,  resonant  modes  of  a 
structure  with  relatively  large  vibratory  displacements  (which  usually  occur  at  lower  order 
modes).  However,  particle  dampers  may  provide  some  damping  under  extreme  conditions 
where  other  options  are  not  feasible. 

The  next  step  in  the  methodology  is  to  select  a  specific  particle  damping  configuration.  To  aid 
in  the  selection  of  an  effective  particle  damper  configuration,  a  list  of  basic  “rules-of-thumb” 
have  been  established.  These  rules  are  based  on  measured  results  from  experimental  testing, 
predicted  results  from  analytical  modeling,  and  published  results  from  other  researchers  [37- 
38].  Key  design  guidelines  include  the  following: 
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•  Total  particle  mass  appears  to  have  a  fairly  significant  effect  on  damping  for  both  single 
particle  and  multiple-particle  dampers.  Increasing  the  mass  tends  to  increase  the 
damping.  Conversely,  decreasing  the  mass  tends  to  reduce  damping.  As  might  be 
expected,  changes  in  the  total  particle  mass  can  lead  to  a  fairly  significant  shift  in  the 
frequency  of  peak  response. 

•  Increases  in  the  coefficient  of  friction  or  viscoelasticity  of  the  particles  tend  to  increase 
the  damping  for  single-particle  dampers.  For  multiple-particle  dampers,  very  little 
effect,  or  perhaps  even  a  slight  adverse  effect,  on  the  damping  is  observed. 

•  For  single-particle  dampers,  orienting  the  dampers  such  that  gravity  is  in  the  direction 
of  the  vibratory  motion  tends  to  reduce  the  damping.  The  orientation  of  gravity  appears 
to  have  much  less  effect  on  dampers  with  multiple  particles. 

•  Both  single-  and  multiple-particle  dampers  appear  to  exhibit  maximum  damping  at  an 
optimum  excitation  level.  Lesser  damping  is  observed  at  excitations  above  and  below 
this  specific  level.  This  relationship  is  closely  linked  to  the  interior  dimensions  of  the 
cavity,  particularly  the  length  of  the  cavity  in  the  direction  of  the  vibratory  motion. 

•  Single-particle  dampers  appear  to  be  more  sensitive  to  changes  in  the  various  particle 
damper  parameters.  This  result  illustrates  a  potential  advantage  of  multiple-particle 
dampers  in  that  precise  tuning  of  the  damper  parameters  may  be  less  critical. 

•  At  times,  both  the  single-  and  multiple-particle  dampers  may  cause  random,  somewhat 
chaotic,  behavior  of  the  damped  system  (e.g.,  the  system  may  exhibit  several  local 
response  peaks  instead  of  a  single  resonant  peak).  This  behavior  occurs  more  commonly 
in  systems  where  the  particle  damper  is  contributing  significant  damping  and  may  result 
due  to  slight  mistuning  of  the  damper  parameters. 

Once  a  specific  particle  damper  configuration  has  been  selected,  the  characteristics  of  the 
undamped  system  again  must  be  determined  to  incorporate  the  effects  of  any  added  (or 
reduced)  mass  due  to  the  particle  damper  cavity.  Often  it  is  not  practical  to  simulate  the  entire 
physical  system  with  the  particle  damper  attached.  Therefore,  the  specific  mode  of  the  system 
to  be  damped  is  modeled  as  a  SDOF,  sinusoidally  excited,  mass-spring-damper  system.  It  is 
necessary  to  determine  the  generalized  mass,  spring  stiffness,  viscous  damping  coefficient,  and 
drive  force  from  the  undamped  physical  system  to  develop  the  analytical  SDOF  model.  Even 
for  cases  where  the  entire  structure  is  to  be  modeled,  it  is  necessary  to  perform  this  step  to 
ensure  that  the  added  cavity  mass  does  not  drastically  affect  the  natural  frequencies  or  modes  of 
interest. 

Once  a  particle  damper  configuration  has  been  selected  and  the  characteristics  of  the  undamped 
system  with  adjusted  mass  have  been  determined,  the  analytical  model  can  be  used  to  evaluate 
the  damping  effectiveness.  Effectiveness  can  be  evaluated  by  comparing  vibratory 
displacements  for  the  undamped  and  damped  cases.  Similarly,  the  relative  effectiveness  of 
particle  damper  designs  can  be  determined  by  comparing  predicted  displacements.  Specific 
details  on  using  the  analytical  particle  damper  model  for  specific  applications  are  given  in 
Sections  8  and  9.  In  these  sections,  predicted  results  from  analytical  particle  damper  models  are 
compared  with  experimental  test  results  for  a  cantilevered  beam  and  an  aluminum  box 
structure. 
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Based  on  the  results  of  the  particle  damper  simulation,  it  may  be  determined  that  the  chosen 
damper  configuration  does  not  provide  sufficient  damping.  Even  if  the  damping  is  sufficient,  it 
may  be  desirable  to  investigate  alternative  designs  to  determine  if  similar  damping  levels  can 
be  achieved  with  lower  added  mass  or  using  more  readily  available  materials.  The  guidelines 
given  above  provide  some  guidance  in  identifying  new  damper  configurations.  If  the  damper 
cavity  configurations  remain  the  same,  new  simulations  can  be  performed  using  the  previously 
determined  system  parameters.  Otherwise,  the  characteristics  of  the  undamped  system  with  the 
new  cavity  configuration  must  be  determined.  Several  design  cycles  may  be  necessary  to  define 
an  effective  particle  damper  configuration. 

Lastly,  once  an  acceptable  particle  damper  configuration  has  been  established,  the  damper 
should  be  fabricated  and  applied  to  the  structure.  Experimental  measurements  should  be  taken 
to  verify  the  effectiveness  of  the  particle  damper  and  to  ensure  that  there  are  no  adverse  effects 
on  other  modes  of  the  system  or  clearance  issues. 
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8  Correlation  with  Experimentai  Beam  Tests 

Experimental  testing  and  analytical  simulations  were  performed  to  validate  the  model  discussed 
in  the  previous  section.  A  cantilevered  beam  was  tested  undamped  and  with  various  impact 
dampers  and  particle  dampers  attached.  The  test  data  was  then  compared  to  analytical 
prediction. 

8.1  Experimental  Testing 

The  first  flexural  mode  of  a  cantilevered  aluminum  beam  was  used  for  the  experimental  testing. 
The  behavior  of  a  vibrating  beam  is  well  understood  and  permits  the  use  of  simple  hand 
calculations  or  finite  element  models  to  evaluate  the  effect  of  various  modifications,  such  as  the 
shift  in  natural  frequency  due  to  added  mass.  This  system  also  can  be  used  to  roughly 
approximate  a  single  degree-of-freedom  system.  The  testing  was  performed  using  an  aluminum 
beam  12.0  inches  in  length,  1.50  inches  wide,  and  0.175  inch  thick,  as  shown  in  Figure  42.  The 
beam  was  rigidly  fixed  at  the  root  and  had  a  first  flexural  mode  of  slightly  less  than  40  hertz. 


Figure  42.  Aluminum  Beam  Used  for  Correlation  Testing 


Experimental  testing  was  performed  on  the  undamped  beam,  on  beams  with  added  mass,  and 
on  beams  damped  with  impact  and  particle  dampers  containing  stainless  steel  spheres.  The 
impact  and  particle  damper  cavities  were  constructed  from  square  acrylic  tubing  with  an  outer 
edge  length  of  0.500  inch  and  a  wall  thickness  of  approximately  0.065  inch.  Aluminum  end 
caps  were  epoxied  to  the  tubing  to  enclose  the  cavities.  Stainless  steel  spheres  of  various  radii 
were  used  for  the  particles. 

8.2  Analytical  Simulation 

Analytically,  the  beam  used  for  the  experimental  testing  is  modeled  as  a  simple  mass-spring- 
dashpot  system,  as  shown  in  Figure  43.  The  beam  is  modeled  using  a  lumped  mass  (at  the 
master  cavity  node)  attached  to  a  damped  spring-to-ground  element.  The  equivalent  mass, 
spring  stiffness,  and  viscous  damping  coefficient  are  chosen  to  simulate  the  undamped  beam. 
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Figure  43.  X3D  Model  Used  to  Simulate  Particle  Dampers 

The  cantilevered  aluminum  beam  used  for  the  experimental  testing  is  12.0  inches  in  length, 
1.50  inches  wide,  and  0.175  inch  thick,  and  has  an  elastic  modulus  of  10.0x10®  psi  and  a 
density  of  0.100  Ibf/in^  (2.59x10'^  Ibf-sVin"^).  The  mass  of  the  beam  can  be  calculated  as 
follows: 

1 W  -f  ^ 

m  =  pV  =  (0.100Xl2.0Xl.50X0.175)=  0.315  Ibf  =  8.15851x10'"  — (34) 


The  equivalent  mass,  ,  for  a  cantilever  beam  of  mass,  m,  carrying  an  end  mass,  M ,  is  given 
as  follows  [39]: 

meq=M  +  0.23m  (35) 


For  the  undamped  beam,  the  equivalent  mass  can  be  found  as  follows: 

-ii  2 

m,,  =  0.23  m  =  (0.23X8. 1585  lx  10'")=  1.87646x10'" 

The  stiffness  of  the  beam,  k ,  is  calculated  as  follows: 


k 


3E1  (3)(l0.0xl0‘)(6.9922xl0-7_^^^,„.,',,„,lM 

1^  (12X  in 


where  the  moment  of  inertia,  I ,  is  calculated  as  follows: 

I  =  —  =  =  6.69922x10'"  in" 

12  12 


The  natural  frequency  of  the  beam,  Wbeam  or  fbeam»  is  calculated  as  follows: 


(36) 


(37) 


(38) 
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°  or  = 39.623  H.  (39) 

The  measured  natural  frequency  of  the  undamped  beam  is  39.37  Hz,  or  247.369  rad/s.  The 
analytical  simulations  use  an  adjusted  stiffness  of  1.14823x10*  Ibf/in  such  that  the  natural 
frequency  is  identical  to  the  measured  frequency  of  39.37  Hz. 

The  beam  system  is  excited  by  a  prescribed  sinusoidal  force  applied  at  the  master  cavity  node. 
Initial  estimates  for  the  analytical  drive  force  are  made  based  on  the  tip  force  equivalent  to  the 
measured  force  input  near  the  root  of  the  beam.  Initial  damping  estimates  are  made  based  on 
the  measured  loss  factor  from  experimental  testing.  The  excitation  force  and  damping  are 
adjusted  such  that  the  predicted  response  of  the  undamped  beam  from  the  analytical 
simulations  matches  the  undamped  response  measured  experimentally.  Thus,  the  analytical 
model  is  tuned  for  the  undamped  response.  Subsequent  damped  analyses  utilize  the  tuned 
model  parameters.  The  drive  force,  equivalent  mass,  spring  stiffness,  and  viscous  damping 
coefficient  used  to  model  the  undamped  beam  in  X3D  for  an  excitation  of  200  mV  RMS  to  the 
shaker  are  given  in  Table  4. 

Table  4.  Parameters  Used  to  Simulate  Undamped  Beam  for  X3D  Analyses 


Parameter 

Value 

Drive  force 

4.84040x10'^  Ibf 

Equivalent  mass 

1.87646x10'^  Ibf-s^/in 

Viscous  damping  coefficient 

3.08293x10'^*  Ibf-s/in 

Spring  stiffness 

1.14823x10*  Ibf/in 

For  damped  simulations,  the  cavity  is  modeled  using  contact  surfaces  defined  by  nodes  which 
are  linked  to  the  master  cavity  node.  Particles  are  tracked  using  a  node  at  the  center  of  each 
particle.  The  particles  are  initially  given  a  random  distribution  within  the  cavity  and  include  a 
gravity  load  which  ensures  that  the  particles  pack  on  the  bottom  of  the  cavity.  The  mass  of  the 
particles  is  included  based  on  the  prescribed  number  of  particles  and  the  radius  and  density  of 
each  particle.  Any  additional  mass  due  to  the  cavity  must  be  added  to  the  equivalent  mass  of 
the  undamped  beam.  Material  properties  are  given  for  the  particles  to  incorporate  the  elastic  or 
viscoelastic  impact  behavior.  A  coefficient  of  friction  also  is  given  for  the  particles  to  account 
for  friction  interactions  between  the  individual  particles  or  between  the  particles  and  the  cavity 
walls. 

A  sample  input  deck  for  the  aluminum  beam  with  a  particle  damper  cavity  with  interior 
dimensions  of  0.372  inch  wide,  0.372  inch  high,  and  0.263  inch  long,  containing  64  0.03125- 
inch  radius  particles  is  given  in  Appendix  A.  The  basic  geometry  of  the  system  is  defined  using 
the  NODE  COORDINATES,  ELEMENTS,  and  CONTACT  cards.  The  position  of  various 
nodes  within  the  model  are  defined  using  the  NODE  COORDINATES  card.  These  nodes  are 
required  to  define  the  cavity  geometry  and  master  cavity  node,  as  well  as  to  track  the  particle 
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centers.  The  damped  spring-to-ground  element  is  created  using  the  ELEMENTS  card  and 
provides  a  link  between  the  master  cavity  node  and  a  fixed  point  (ground).  The  spring  stiffness 
and  viscous  damping  coefficient  for  this  element  are  defined  on  the  MAT  ID  card.  The 
CONTACT  card  is  used  to  define  12  triangular  master  contact  surfaces,  which  define  the  cavity 
geometry,  along  with  slave  nodes  corresponding  to  the  center  of  each  particle.  Properties  for 
the  particles  are  defined  on  the  PDMP  card.  Basic  input  parameters  for  the  PDMP  card  are 
given  in  Appendix  B  (the  PDMP  card  has  been  created  specifically  for  particle  damper 
simulations  and  is  not  included  in  the  basic  X3D  user’s  manual).  The  BODY  card  is  used  to 
simulate  gravity  such  that  the  particles  tend  to  fall  to  the  bottom  of  the  cavity.  The  equivalent 
mass  of  the  system  is  lumped  at  the  master  cavity  node  using  the  MASS  card.  Note  that  an 
equivalent  mass  of  2.04213x10  '*  Ibf-s^/in  is  used  which  includes  the  mass  of  the  aluminum 
beam  plus  the  additional  mass  of  the  damper  cavity  (but  not  the  particles).  The  RBEL  card 
links  the  cavity  nodes  to  the  master  cavity  node  and  is  used  to  ensure  that  the  entire  cavity 
moves  as  a  single  unit.  The  BOUNDARY  CONDITIONS  card  ensures  that  the  cavity  only 
moves  in  the  direction  of  the  applied  force.  The  TITLE  card  allows  a  descriptive  title  to  be 
assigned  to  the  analysis.  The  PARAMETERS  card  sets  parameters  relating  to  the  total 
simulated  time,  the  timestep,  the  total  number  of  timesteps,  and  the  frequency  at  which  results 
are  output  to  the  trace  file.  For  these  analyses,  simulations  are  run  for  15  seconds  using  a 
timestep  of  0.5  microsecond.  Results  are  output  to  the  trace  file  every  millisecond.  Lastly,  the 
TRACE  card  is  used  to  define  those  nodes  for  which  output  to  the  trace  file  is  requested. 
Typically,  results  for  the  master  cavity  node  are  output  along  with  results  for  each  of  the 
particles.  Further  details  on  the  various  options  available  for  each  of  these  cards  (with  the 
exception  of  the  PDMP  card)  are  given  in  the  X3D  user’s  manual. 

The  excitation  force  used  for  a  simulation  is  defined  in  the  USRFRC  subroutine  which  must  be 
included  in  the  compilation  of  the  X3D  code.  The  X3D  code  contains  a  skeleton  version  of  the 
USRFRC  routine.  A  sample  USRFRC  subroutine  modified  for  particle  damper  simulations  is 
given  in  Appendix  C.  This  sample  routine  imposes  a  force  of  4.84040x10'^  Ibf  at  a  frequency  of 
37.21  Hz,  or  233.797  rad/s.  The  force  is  applied  at  node  25,  the  master  cavity  node.  To  perform 
simulations  at  other  drive  frequencies  or  excitation  levels,  or  with  other  master  cavity  nodes, 
this  subroutine  must  be  modified  and  a  new  version  of  X3D  compiled.  The  master  cavity  node 
also  is  currently  hard-coded  on  two  lines  (which  set  the  variable  NODEM)  in  the  PDFORC 
routine.  These  lines  also  must  be  modified  to  change  the  master  cavity  node. 

8.3  Results 

Figure  44  shows  selected  frames  from  early  in  the  64  particle  simulation  as  discussed  above. 
Note  that  the  particles  are  initially  given  a  random  distribution  within  the  cavity  (as  shown  in 
the  upper  left  frame  at  a  time  of  0.0  s)  but  quickly  pack  on  the  bottom  of  the  cavity  due  to 
gravity  (as  shown  in  the  lower  right  frame  at  a  time  of  0.10  s).  The  frames  shown  below  have 
been  created  using  a  MATLAB  routine  which  processes  the  data  output  to  the  trace  file.  The 
frames  can  be  animated  to  observe  the  cavity  and  particle  motion. 
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Figure  44.  Selected  Frames  from  Particle  Damper  Simulation 

Various  simulations  have  been  performed  corresponding  to  cases  which  have  been  tested 
experimentally.  Measurements  have  been  made  for  the  undamped  beam,  the  beam  with  an 
added  mass  equal  to  the  mass  of  the  dampers,  and  for  various  particle  damper  configurations 
(of  the  same  total  particle  mass)  containing  1  0.250-inch  diameter  particle,  64  0.0625-inch 
diameter  particles,  and  512  0.03125-inch  diameter  particles.  Figure  45  shows  experimental  and 
analytical  beam  tip  displacements  at  various  frequencies  for  the  undamped,  added  mass,  and 
single  particle  configurations  under  a  force  equivalent  to  a  200  mV  RMS  excitation  signal  to 
the  shaker.  Figure  46  shows  results  for  the  various  damped  cases  under  a  force  equivalent  to  a 
200  mV  RMS  excitation  signal  to  the  shaker.  Similar  results  corresponding  to  a  400  mV  RMS 
excitation  signal  to  the  shaker  are  shown  in  Figure  47  and  Figure  48. 
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For  all  of  the  damper  analyses,  the  particles  are  given  viscoelastic  material  properties  with  Eq 
equal  to  27.5x10®  psi,  Ei  equal  to  5.0x10®  psi,  and  Ti  equal  to  2.0x10'®  second  (corresponding 
to  a  maximum  equivalent  loss  factor  of  0.083).  A  coefficient  of  friction  of  0.30  has  been  used 
in  the  analyses.  For  some  of  the  experimental  results,  two  sets  of  data  are  given  for  the  same 
configuration.  In  each  case,  the  first  set  of  results  corresponds  to  data  taken  with  the  excitation 
frequency  increasing  during  testing  and  the  second  set  with  the  excitation  frequency 
decreasing.  An  interesting  result  is  observed  in  the  experimental  data  for  the  damper  with  512 
particles.  Under  200  mV  RMS  excitation,  very  little  attenuation  is  seen.  However,  under  400 
mV  RMS  excitation,  considerable  attenuation  is  observed  once  the  displacements  reach  a 
certain  level  (whether  the  excitation  frequency  is  increasing  or  decreasing).  During 
experimental  testing  it  has  been  observed  that  the  small  particles  in  this  damper  tend  to  clump 
together,  and  it  is  believed  that  there  may  be  cohesive  forces  between  the  particles  due  to  static 
electric  charges  or  oil  on  the  particles  which  influence  the  results. 

In  general,  the  particle  damper  simulation  code  predicts  the  correct  trends.  For  example,  the 
model  predicts  little  attenuation  for  the  single  particle  dampers,  but  considerably  more 
attenuation  for  multiple-particle  dampers  containing  64  and  512  particles.  The  model  also 
captures  the  experimental  trend  of  optimum  damping  at  a  specific  number  of  particles  with 
lesser  attenuation  when  either  a  smaller  or  greater  number  of  particles  is  used.  This  trend  is 
evident  in  the  increased  attenuation  observed  when  the  number  of  particles  is  increased  from  a 
single  particle  to  64  particles;  and  the  decreased  attenuation  observed  when  the  number  of 
particles  is  further  increased  to  512  particles.  The  model  does  tend  to  overpredict  the 
attenuation,  particularly  for  multiple-particle  dampers. 
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Figure  45.  Beam  Tip  Displacements  with  200  mV  RMS  Excitation  Force 
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Figure  46.  Damped  Tip  Displacements  with  200  mV  RMS  Excitation  Force 
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Figure  47.  Beam  tip  displacements  with  400  mV  RMS  excitation  force 
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Figure  48.  Damped  Tip  Displacements  with  400  mV  RMS  Excitation  Force 


9  Correlation  with  Experimental  Chassis  Tests 

Additional  experimental  testing  and  analytical  simulations  have  been  performed  on  an 
aluminum  chassis.  The  basic  experimental  test  setup,  analytical  modeling  details,  and 
correlation  of  experimental  results  with  analytical  predictions  are  discussed  in  the  following 
paragraphs. 

9.1  Experimental  Testing 

The  aluminum  chassis  used  for  this  testing,  along  with  corresponding  test  hardware,  is  shown 
in  Figure  49.  The  structure  and  related  hardware  are  from  a  previous  effort  and  provide  a 
suitable  article  for  particle  damping.  The  chassis  is  basically  an  aluminum  box  with  dimensions 
of  17.0  inches  in  length,  8.0  inches  in  width,  and  3.0  inches  in  depth.  The  box  is  constructed 
from  aluminum  sheet  with  a  nominal  thickness  of  0.049  inch.  The  box  is  supported  by  a 
relatively  stiff  aluminum  structure  as  shown  in  the  figure.  The  entire  lower  edge  of  the  box  is 
fixed. 


Figure  49.  Aluminum  Chassis  Structure  and  Test  Hardware 

The  particle  dampers  are  focused  on  attenuating  the  fundamental  mode  of  the  structure.  For  this 
mode,  most  of  the  motion  is  seen  in  the  upper  surface  of  the  box  with  the  displacements 
following  roughly  half  of  a  sinusoidal  cycle  in  both  the  width  and  length  directions.  The 
deformed  shape  of  the  first  fundamental  mode,  as  predicted  by  finite  element  analysis,  is  shown 
in  Figure  50. 
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Figure  50.  Fundamental  Mode  of  Box  Structure 

Similar  to  the  beam  testing,  the  box  is  excited  through  a  stinger  which  is  centered  across  the  8.0 
inch  width  of  the  panel  and  located  3.425  inches  above  the  bottom  of  the  panel.  A  force  gage  is 
located  between  the  stinger  and  the  box  to  record  the  actual  force  input  into  the  structure.  A 
sinusoidal  signal  is  passed  to  the  shaker,  with  the  input  voltage  controlled  to  maintain  a  specific 
excitation  force.  An  accelerometer  is  positioned  at  the  center  of  the  panel,  on  the  inner  surface 
of  the  box,  to  record  accelerations. 

The  damper  is  placed  at  the  center  of  the  panel,  on  the  outer  surface  of  the  box,  where  the 
largest  displacements  are  seen  for  this  mode  (and  are  relatively  large  compared  to  the 
maximum  displacements  of  other  higher  order  modes).  The  damper  cavity  is  fabricated  from 
square  brass  tubing  with  an  outside  edge  length  of  0.500  inch  and  a  wall  thickness  of  0.028 
inch,  yielding  an  interior  edge  length  of  0.444  inch.  End  caps  constructed  from  a  0.028-inch- 
thick  brass  sheet  are  epoxied  to  the  ends  of  the  tubing  to  create  the  cavities.  The  cavities  have 
an  interior  length  of  0.350  inch.  Initial  studies  have  been  performed  using  particle  dampers 
containing  128  spherical  stainless  steel  particles  with  a  diameter  of  0.0625  inch. 

9.2  Analytical  Simulation 

Analytical  particle  damper  simulations  have  been  performed  using  both  a  SDOF  model  and  a 
multiple  degree-of-freedom  (MDOF)  model.  The  SDOF  model  is  similar  to  the  model  used  for 
the  beam  simulations.  The  aluminum  box  structure  is  modeled  as  a  simple  mass-spring-dashpot 
system,  with  the  analytical  model  tuned  to  match  the  response  of  the  undamped  structure.  The 
drive  force,  equivalent  mass,  spring  stiffness,  and  viscous  damping  coefficient  used  to  model 
the  undamped  box  structure  in  X3D  for  a  0.1  Ibf  RMS  excitation  are  given  in  Table  5. 

Table  5.  Parameters  used  to  simulate  undamped  box  structure  for  X3D  analyses 


Parameter 

Value 

Drive  force 

5.65771x10'^  Ibf 
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Equivalent  mass 

3.54282x10'^  Ibf-s^/in 

Viscous  damping  coefficient 

3.55000x10'^  Ibf-s/in 

Spring  stiffness 

1.95442x10^  Ibf/in 

. 

Additional  analytical  simulations  have  been  performed  using  a  MDOF  model.  For  these 
simulations,  the  entire  box  structure  is  modeled.  A  relatively  eoarse  mesh  is  used  since  only  the 
fundamental  mode  is  of  concern  and  to  reduce  the  likelihood  of  higher  frequency  noise  in  the 
solution.  The  model  used  for  these  analyses  is  shown  in  Figure  51.  For  these  simulations,  the 
force  from  the  stinger  is  applied  at  the  actual  location  used  for  the  experimental  testing,  rather 
than  at  the  master  cavity  node.  The  node  at  the  center  of  the  panel  is  used  as  the  master  cavity 
node  and  the  predicted  vibratory  accelerations  are  recorded  at  this  node.  The  cavity  and 
particles  are  modeled  and  tracked  as  in  the  SDOF  simulations. 


Figure  51.  Model  Used  for  MDOF  Simulations 


9.3  Results 

Experimental  testing  and  analytical  simulations  have  been  performed  at  various  excitation 
levels  and  across  a  range  of  frequencies.  The  response  of  the  undamped  chassis  structure,  and 
the  ehassis  structure  with  particle  dampers  eontaining  128  particles,  has  been  investigated. 
Results  are  presented  as  accelerance  (in  this  case,  the  ratio  of  the  vibratory  accelerations  at  the 
center  of  the  panel  to  the  input  force  at  the  shaker  location)  versus  frequency.  Figure  52 
compares  results  from  the  experimental  testing  and  simulations  using  the  SDOF  model.  Some 
interesting  results  are  seen  in  the  experimental  data.  At  the  0.1  Ibf  RMS  excitation  level,  the 
peak  response  of  the  damped  chassis  structure  has  approximately  the  same  magnitude  as  the 
undamped  structure,  but  occurs  at  a  slightly  higher  frequency.  As  the  excitation  level  is 
increased  up  to  0.5  Ibf  RMS,  the  peak  response  is  attenuated  considerably  with  the  frequency  of 
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the  peak  increasing  with  increasing  excitation.  Above  0.5  Ibf  RMS  excitation  the  damping 
decreases,  but  the  frequency  of  the  peak  continues  to  increase.  It  is  believed  that  this  behavior 
results  due  to  a  combination  of  various  effects,  including  changes  in  damping  and  effective 
added  mass.  Such  behavior  is  not  seen  in  the  analytical  results.  The  analyses  predict  increasing 
damping  with  increasing  excitation,  with  the  frequency  of  the  peak  response  decreasing  with 
increasing  excitation. 
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Figure  52.  Comparison  of  Experimental  Measurements  and  SDOF  Analysis  Results 

Analyses  with  the  MDOF  model  have  been  performed  to  determine  if  the  discrepancy  between 
the  experimental  and  SDOF  model  results  is  due  to  the  inability  of  the  SDOF  model  to 
accurately  capture  the  dynamics  of  the  chassis  structure.  Figure  53  compares  results  from  the 
experimental  testing  and  simulations  using  the  MDOF  model.  Results  from  the  MDOF  model 
are  similar  to  those  from  the  SDOF  model.  Although  the  undamped  model  could  be  more 
accurately  tuned  to  match  the  experimental  response  of  the  undamped  structure,  it  is  unlikely 
that  this  would  have  much  of  an  effect  on  the  results.  Changes  in  the  coefficient  of  friction  and 
viscoelastic  properties  do  not  appear  to  have  a  large  effect  on  the  overall  peak  response.  As  a 
result,  it  is  thought  that  the  momentum  transfer  may  have  a  critical  effect  on  the  damping. 
Reductions  in  the  forces  transferred  to  the  structure  reduce  the  damping  and  correlate  better 
with  the  experimental  results.  It  is  possible  that  compliance  of  the  cavity  walls,  which  are 
currently  assumed  to  be  rigid,  may  contribute  to  the  discrepancy.  It  is  recommended  that  future 
efforts  focus  on  ensuring  that  appropriate  forces  are  applied  during  impacts  between  the 
particles  and  the  cavity  walls. 
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Figure  53.  Comparison  of  Experimental  Measurements  and  MDOF  Analysis  Results 
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10  Proof  of  Concept 

Engine  hardware  for  application  of  the  MPID  treatment  using  the  design  methodology  was 
selected  with  the  aid  of  Pratt  &  Whitney.  The  hardware  was  a  divergent  flap  backstructure 
subcomponent  for  the  high  speed  civil  transport.  The  four-ribbed  subcomponent  shown  in 
Figure  54  is  approximately  12  by  12  by  4  inches  high.  A  two-ribbed  structure  has  been  made 
available  to  CSA  for  design  and  testing  of  particle  dampers.  The  structures  were  designed  for 
high  temperature  capability  and  acoustic  noise  reduction.  They  were  made  from  gamma- 
titanium  to  reduce  part  weight.  This  material  is  also  brittle  and  has  very  low  inherent  damping. 
High  cycle  fatigue  is  a  large  problem  for  these  types  of  structures. 


Figure  54.  Four-Ribbed  Subcomponent  at  Pratt  &  Whitney  Test  Facility 

A  finite  element  model  of  the  two-ribbed  structure  was  created.  Both  ribs  were  modeled  using 
finite  element  (FE)  and  are  shown  in  Figure  55.  The  first  four  distinct  modes  are  shown  in 
Figure  56. 
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was  used  with  a  1,000-Hz  bandwidth.  The  first  mode  is  apparent  at  177  Hz.  The  second  mode 
is  240  Hz.  The  first  mode  is  the  mode  to  damp.  Particle  damping  will  have  effect  on  the  higher 
modes,  as  well. 

The  first-mode  peak  displacement  of  the  undamped  structure  was  3.25e-4  inches.  An  impact 
damper  with  slightly  less  than  6.5e-4  clearance  between  the  particle  and  cavity  wall  should 
have  the  most  influence  on  the  structure. 


r'"' 


Figure  57.  Substructure  Mounted  on  the  Shaker  Table  (1  50-ParticIe  Damper  Attached) 
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Figure  58.  FRF  for  Undamped  Structure 

The  FE  model  was  used  to  calculate  the  moments  of  inertia.  As  a  representative  aircraft 
structure,  weight  was  considered  important.  A  simulation  was  run  for  a  MPID  containing  50 
0.0625-inch-diameter  stainless  steel  302  particles.  The  cavity  was  brass  having  a  0.25-inch 
inner-length  cubic  interior  with  0.027-inch  wall  thickness.  A  comparison  of  the  simulated 
results  with  the  undamped  structure  is  shown  in  Figure  59.  This  run  showed  that  the  MPID 
would  give  good  results,  but  it  over-predicted  the  damping.  The  measured  FRF  for  the  single 
cavity  with  50  particles  is  shown  in  Figure  60. 
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To  damp  the  outer  portions  of  the  structure,  three  0.25-inch  inner-length  cubic  cavities  were 
constructed  and  each  was  filled  with  50  stainless  steel  302  alloy  spherical  particles  having 
0.0625-inch  diameter.  The  cavities  were  applied  at  the  center  and  outer  “wings,”  as  shown  in 
Figure  61. 


Figure  61.  Substructure  with  3  MPIDs  Attached 

The  FRF  for  three  cavities  is  shown  in  Figure  62.  As  expected,  the  application  of  the  three 
cavities  affected  the  both  the  mode  at  177  Hz  and  the  mode  at  240  Hz. 

The  simulation  code  was  used  to  determine  that  given  the  material  and  geometric  constraints  of 
cavity  material  and  size,  particle  material  and  size,  and  structure  dynamics,  102  particles  would 
give  better  damping  than  50  particles.  Three  cavities  with  102  0.0625-diameter  particles  were 
constructed  and  applied.  The  FRF  is  shown  in  Figure  63.  The  increase  in  the  number  of 
particles  had  a  relatively  large  effect  on  both  the  first  and  second  modes.  The  first  mode  was 
reduced  by  23  dB.  The  second  mode  as  measured  by  the  right-hand  accelerometer  was  reduced 
by  49  dB. 
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Figure  62. 3  50-Particle  Dampers 


Tliree  102-particle  dampers,  2.5V  max  load 


Figure  63. 3  102-Particle  Dampers 


11  Commercialization 

Commercialization  is  an  important  part  of  the  SBIR  program.  CSA  applied  many  of  the  lessons 
learned  during  this  SBIR  to  design  a  particle  damping  treatment  for  a  commercial  project. 

11.1  Introduction 

A  commercial  jet-engine  manufacturer  is  developing  and  integrating  an  engine  that  on  occasion 
experiences  harmful  vibration.  This  vibration  is  induced  upon  stall,  a  rare  event.  The  level  of 
vibration  is  unacceptable.  Among  other  things,  it  causes  excessive  wear.  The  vibration  is 
predominantly  in  a  mode  with  a  frequency  of  about  120  Hz  in  the  test  stand  configuration.  This 
frequency  is  expected  to  increase  to  roughly  144  Hz  in  the  shipboard  application  with  a 
different  mounting  configuration  and  various  test  components  removed.  The  engine  casing 
operates  at  a  relatively  high  temperature,  in  the  range  of  600-800  T.  Vibration  damping  using 
ordinary  methods  is  not  feasible  at  this  temperature. 

The  customer  would  like  to  eliminate  the  vibration  problem.  One  potentially  effective  approach 
is  to  add  damping  to  the  system.  CSA  Engineering  is  a  company  that  specializes  in  vibration 
suppression.  CSA  has  several  available  methods  for  damping  vibration,  although  many  of  these 
are  not  applicable  at  elevated  temperature. 

CSA  proposed  to  employ  alternative  methods  for  reducing  the  harmful  vibration  to  a  point  at 
which  its  level  is  acceptable.  In  this  development  effort,  CSA  investigated  MPID.  CSA 
designed  concepts  and  produced  a  prototype  for  inclusion  in  upcoming  testing. 

11.2  Design 

The  following  data  were  supplied  by  the  customer: 

1.  Typical  case  vibration  is  roughly  0.4  inch  per  second  (ips)  rms.  Max  case  vibration 
during  stall  condition  reaches  3.0  ips  rms. 

2.  The  entire  case  weighs  1950  lb.  The  modal  mass  participation  in  the  146  Hz  mode  of 
concern  is  600  lb.  The  mode  shape  is  shown  in  Figure  64. 

3.  During  the  stall,  the  vibration  increased  from  0.4  ips  rms  to  3.0  ips  rms  in  slightly  less 
than  5  seconds  (Figure  65). 

4.  The  measured  damping  ratio  for  the  mode  is  0.067. 
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Figure  64.  Engine  Mode  Shape 


— —  visa  X'«13:ZS:D2.iI3  ¥<1.2932  -  VMCV  [ISS-IWSI  30<i3:2fi:Di.it3  TC«3.i6313 

—  vaar  igg-ws) 


Figure  65.  Vibration  at  Stall 

The  peak  displacement  was  calculated  as  0.0046  (zero-peak)  inche  for  the  given  modal  mass, 
damping,  stiffness,  and  velocity.  The  relatively  high  damping  ratio  of  0.067  required  that  the 
particle  damper  be  tuned  to  produce  maximum  damping  for  the  modal  mass,  frequency,  and 
displacement  of  the  mode. 

To  obtain  maximum  damping,  the  particles  in  the  particle  damper  must  impact  the  ends  of  the 
cavity.  With  the  given  displacement,  the  cavity  free  length  had  to  be  approximately  0.009 
inche.  Also,  for  maximum  effect,  the  particles  should  not  interfere  with  each  other.  This  is 
accomplished  with  separate  cavities  for  each  particle.  This  is  known  as  a  multiple  individual 
particle  (MIP)  configuration.  It  is  essentially  many  side-by-side  impact  dampers. 

Stainless  steel  type  302  (SS-302)  was  chosen  for  the  particle  and  cavity  material.  CSA  has 
found  that  SS-302  provides  higher  damping  than  most  other  particle  materials,  withstands  high 
temperature,  and  shows  little  wear  after  long  duration. 

Two  reconfigurable  particle  dampers  were  designed  that  had  separate  cavities  for  each  particle. 
Shims  were  used  to  adjust  the  length  of  the  cavities.  The  dampers  consisted  of  four  10  by  10- 
inch  aluminum  plates  with  thickness  equal  to  the  particle  diameters.  Aluminum  was  used  in 
prototyping,  but  for  high-temperature  application,  the  plates  should  be  manufactured  from  the 
same  material  as  the  particles,  i.e.,  stainless  steel.  The  four  plates  were  drilled  to  hold  particles 
and  were  sandwiched  between  0.050-inch-thick  end  plates.  The  general  plate  configuration  is 
shown  in  Figure  66.  Shims  were  placed  between  the  damper  plates  and  end  plates  to  vary  the 
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cavity  free  lengths.  Two  different  particle  sizes  were  examined:  0.25-  and  0.5-inch  diameter. 
The  cavity  plates  for  these  different  particle  sizes  are  shown  in  Figure  67  and  Figure  68.  There 
are  988  particles  in  each  0.025-inch  plate  and  221  particles  in  each  0.5-inch  plate. 


Figure  66.  Particle  Damper  Configuration 
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Figure  67. 0.25-inch  Particle  Cavity  Plate 
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Figure  68. 0.5-inch  Particle  Cavity  Plate 

11.3  Tests 

A  test  apparatus  was  designed  to  approximate  the  modal  properties  of  the  turbine  engine.  For 
safety  and  manageability,  a  modal  mass  of  350  pounds  was  used.  As  such,  the  setup  represents 
a  scaled  model  of  the  mode  and  its  600  pounds  of  modal  mass.  The  test  apparatus  is  shown  in 
Figure  69  and  Figure  70.  Measurements  were  taken  with  an  MTS  661.19E-04  5500-lbf  load 
cell  and  three  Kistler  8730A500M1  500-g  accelerometers.  The  accelerometers  were  evenly 
distributed  around  the  lower  perimeter  of  the  main  mass.  The  output  of  the  three  accelerometers 
was  then  averaged  to  obtain  the  acceleration  of  the  total  mass  in  the  vertical  direction.  Their 
individual  output  was  used  to  verify  that  the  mode  shape  was  the  vertical  plunge  mode. 
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Figure  69.  Test  Apparatus  Model 


Figure  70.  Test  Apparatus  with  0.5-inch  Particle  Cavity  Plate 


The  first  mode  of  the  test  apparatus  at  low  levels  of  excitation  was  132  Hz.  The  damping  ratio 
was  0.047.  The  stiffness  of  the  system  was  6.25e5  Ibf/in.  It  was  calculated  that  to  get  3  ips  rms, 
given  the  measured  damping,  a  force  of  300  Ibf  would  be  required.  It  turned  out,  however,  that 
the  inherent  system  damping  was  amplitude  dependent.  At  1000  Ibf,  the  damping  ratio  was 
0.071  and  the  modal  frequency  was  127  Hz.  At  this  level  of  damping,  it  would  take  a  force  of 
450  Ibf  to  achieve  the  same  velocity.  Also,  the  test  apparatus  was  not  a  true  SDOF  system. 
Therefore,  more  force  was  necessary.  A  series  of  baseline  measurements  were  made.  Swept- 
sine  tests  were  performed  in  which  the  force  amplitude  was  kept  constant  and  the  frequency 
was  incrementally  increased.  The  test  was  then  repeated  with  increased  force  amplitude.  Table 
6  shows  how  baseline  damping  and  modal  frequency  changed  with  amplitude.  Response  plots 
from  the  swept-sine  tests  are  in  the  Appendix. 

Table  6.  Modal  Properties  of  the  Baseline  System  at  Given  Excitation  Levels 


Force  (Ibf) 

Frequency  (Hz) 

Zeta 

150 

132 

0.047 

300 

131 

0.057 

400 

130 

0.060 

500 

130 

0.065 

79 


600 

129 

0.071 

800 

128 

0.071 

1000 

127 

0.071 

1200 

127 

0.070 

There  are  several  mechanisms  for  particle  damping.  They  are  momentum  transfer,  elastic 
restitution  of  the  particles  and  the  cavity,  and  friction.  Multiple  particles  in  a  single  cavity  may 
have  influence  across  a  broad  frequency  bandwidth,  but  they  typically  have  less  overall 
damping  influence  than  single-particle  impact  dampers.  The  highest  levels  of  damping  are  seen 
when  the  particles  get  out  of  phase  with  the  base  motion.  This  is  called  “turning  on.”  This 
phenomenon  can  be  seen  in  the  right-hand  plot  in  Figure  71. 


un<lamped  vs.  O.S  dia balb,  w/ 9mn  shims,  iDrca=299.75  Cbl) 


Figure  71.  Particle  Damping  “Turning  On”  in  the  Graph  on  the  Right 

As  can  be  seen  in  Figure  71,  particle  damping  behavior  can  be  very  sensitive  to  cavity  length. 
To  study  this,  shims  of  0.002,  0.004,  0.006,  0.009,  and  0.012  inch  were  placed  between  the 
cavity  plates  and  the  end  plates.  The  cavity  plates  were  measured  to  be  0.250,  0.251,  0.251, 
0.252, 0.501, 0.5015, 0.502,  and  0.503  inch. 

11.4  Summary 

Calculation  showed  that  0.5-inch  particles  with  0.009-inch  cavity  free  length  would  perform 
best  given  the  turbine  engine  system  parameters. 

Several  observations  were  gained  from  the  test  apparatus: 

1.  The  0.5-inch-diameter  particles  showed  better  performance  than  the  0.25-inch-diameter 
particles  in  all  cases.  This  is  probably  due  to  having  more  moving  mass  than  the  0.25- 
inch  particles:  16.5  lb  versus  9.2  lb. 

2.  In  many  cases,  the  0.25-inch-diameter  particles  provided  some  damping  but  the 
particles  never  “turned  on.” 

3.  The  0.009-inch  shims  gave  slightly  better  performance  than  the  0.006-inch  shims.  Since 
the  0.5-inch  plates  were  0.001  to  0.003  inch  thicker  than  the  0.5-inch  nominal  thickness, 
the  cavity  free  lengths  were  actually  0.010  to  0.012  inch  and  0.007  to  0.009  inch  for  the 
two  cases,  respectively.  The  best  performance  is  likely  somewhere  in-between. 
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When  particles  turn  on,  it  is  practically  impossible  to  use  a  curve  fit  to  derive  a  number  for 
damping.  The  response  is  highly  nonlinear.  Reduction  in  response  is  typically  used  as  the 
measure  of  damping.  A  plot  of  the  reduction  in  velocity  response  where  the  baseline  test  system 
is  closest  in  response  to  the  turbine  engine  is  shown  in  Figure  72.  The  reduction  in  response  for 
the  0.009-inch  shims  is  approximately  35  percent  at  the  peak.  This  was  measured  for  the  test 
apparatus  that  exhibited  higher  inherent  damping  than  the  turbine  engine.  The  reduction  would 
be  higher  for  the  turbine  engine. 

The  impacting  particles  disrupt  the  response  in  this  particular  vibration  mode,  and  potentially  in 
other  modes  that  are  excited.  Because  the  damping  mechanism  is  not  effective  for  extremely 
small  motions,  it  will  be  difficult  to  test  upon  installation.  Fundamentally,  this  damping 
mechanism  acts  as  a  motion  limiter.  A  small  amount  of  vibration  may  occur,  but  the  particular 
damper  parameters  have  been  selected  to  cause  the  damping  to  “turn  on”  and  become  effective 
at  a  relatively  low  vibration  level. 


undamped  vs.  0.5  di a  balls,  w/  9mil  shims,  force=800.25  (Ibf) 


Figure  72.  Reduction  in  peak  velocity  response 

The  particle  dampers  for  the  test  system  used  aluminum  for  the  cavity  material.  For  the 
engine’s  operating  environment,  the  final  system  should  use  SS-302.  Aluminum  was  used  in 
the  test  system  to  keep  machining  costs  down.  There  should  not  be  any  appreciable  change  in 
performance,  however.  The  cavity  wall  thickness  may  be  reduced  to  save  weight,  as  can  the 
surrounding  metal  in  the  cavity  plates.  Standard  practices  for  fatigue  and  strength  should  be 
used. 
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CSA  envisions  two  particle  dampers  placed  on  opposite  sides  of  the  engine  for  maximum 
effect.  The  cavities  should  be  oriented  so  that  their  free  length  is  parallel  to  the  radial  direction 
of  the  engine.  The  overall  size  may  be  altered  for  packaging  and  adaptation  to  the  engine.  If  the 
size  or  number  of  particles  is  altered,  however,  there  will  likely  be  a  commensurate  change  in 
damping. 
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12  Future  Research 

The  work  in  this  report  demonstrated  the  difficulty  of  designing  particle  dampers  for 
engineering  structures.  The  design  capability  is  still  not  easy  to  use  and  is  subject  to  errors. 

Recommendations  for  the  further  study  include  the  following: 

•  Continue  to  attempt  to  resolve  discrepancies  between  code  predictions  and 
experimentally  observed  particle  damping  behavior. 

•  Streamline  code  to  support  analysis  of  MDP  configurations. 

A  strong  basis  in  the  development  of  the  analytical  capability  to  predict  the  performance 
particle  damping  configurations  has  been  established.  Further  work  is  needed  in  terms  of 
correctly  capturing  certain  key  parameters  of  particle  damping  behavior  such  as  the  gradual 
turn  off  of  the  multiple  particle  damper  configurations.  The  model  should  also  be  adapted  to 
work  with  multiple  sizes  of  particles  within  one  cavity. 
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13  Conclusions 

Due  to  the  complex  interactions  of  the  loss  mechanisms  in  a  particle  damper  and  the  large 
number  of  parameters  affecting  the  damper  performance,  it  is  extremely  difficult  to  explicitly 
define  a  particle  damper  configuration  for  a  particular  application.  However,  based  on  the 
damper  behavior  observed  in  experimental  testing  and  analytical  simulations,  design  guidelines 
have  been  established. 

Particle  size  and  number  of  particles  affect  the  response  amplitude  at  which  peak  damping 
occurs.  For  this  particular  cavity  configuration,  the  fewer,  larger  particles  tend  to  peak  at  lower 
amplitudes  (0  to  -5  dB)  while  the  more  numerous,  smaller  particles  favor  the  higher  amplitudes 
(+3  to  +5  dB).  The  damper  with  the  smaller  particles  also  exhibits  a  narrower  peak,  but  with 
higher  peak  damping  values.  In  addition,  the  smaller  particles  are  seen  to  contribute  damping  at 
higher  amplitudes  (+20  to  +30  dB)  whereas  the  larger  particles  yield  almost  no  damping. 
Similar  trends  were  observed  for  a  range  of  particle  materials. 

Based  on  the  experimental  testing  and  analytical  studies,  a  list  of  basic  “rules-of-thumb”  has 
been  established  to  aid  in  the  selection  of  an  appropriate  particle  damper  configuration.  Key 
design  guidelines  include  the  following: 

•  Total  particle  mass  appears  to  have  a  fairly  significant  effect  on  damping  for  both  single 
particle  and  multiple  particle  dampers.  Increasing  the  mass  tends  to  increase  damping. 
Conversely,  decreasing  the  mass  tends  to  reduce  damping.  As  might  be  expected, 
changes  in  the  total  particle  mass  can  lead  to  a  fairly  significant  shift  in  the  frequency  of 
peak  response. 

•  Increases  in  the  coefficient  of  friction  or  viscoelasticity  of  the  particles  tend  to  increase 
the  damping  for  single  particle  dampers.  For  multiple-particle  dampers,  very  little 
effect,  or  perhaps  even  a  slight  adverse  effect,  on  the  damping  is  observed. 

•  For  single  particle  dampers,  orienting  the  dampers  such  that  gravity  is  in  the  direction  of 
the  vibratory  motion  tends  to  reduce  the  damping.  The  orientation  of  gravity  appears  to 
have  much  less  effect  on  dampers  with  multiple  particles. 

•  Both  single-  and  multiple-particle  dampers  appear  to  exhibit  maximum  damping  at  an 
optimum  excitation  level.  Lesser  damping  is  observed  at  excitations  above  and  below 
this  specific  level.  This  relationship  is  closely  linked  to  the  interior  dimensions  of  the 
cavity,  particularly  the  length  of  the  cavity  in  the  direction  of  the  vibratory  motion. 

•  Single-particle  dampers  appear  to  be  more  sensitive  to  changes  in  the  various  particle 
damper  parameters.  This  result  illustrates  a  potential  advantage  of  multiple-particle 
dampers  in  that  precise  tuning  of  the  damper  parameters  may  be  less  critical. 

•  Within  a  large  range  (-94  to  850  deg  F)  particle  damping  is  relatively  temperature 
insensitive.  There  was  a  slight  increase  in  damping  at  higher  temperatures. 

•  At  times,  both  the  single-  and  multiple-particle  dampers  may  cause  random,  somewhat 
chaotic  behavior  of  the  damped  system  (e.g.,  the  system  may  exhibit  several  local 
response  peaks  instead  of  a  single  resonant  peak).  This  behavior  occurs  more  commonly 
in  systems  where  the  particle  damper  is  contributing  significant  damping  and  may  result 
due  to  slight  mistuning  of  the  damper  parameters. 
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The  X3D  explicit  finite  element  code  was  modified  to  perform  particle  damper  simulations. 
The  simulations  are  based  on  the  particle  dynamics  method.  The  utility  of  the  particle  dynamics 
method  is  based  on  the  ability  to  simulate  contact  interaction  forces  between  the  individual 
particles  and  the  cavity  walls,  using  force-displacement  relations  with  a  small  number  of 
parameters  which  capture  the  most  important  contact  properties.  The  following  is  a  summary  of 
the  primary  assumptions  and  limitations  of  the  particle  damper  simulation  code: 

•  All  particles  are  assumed  to  be  identical;  particles  are  assumed  to  be  spherical  in  shape 
with  a  given  radius  and  material  properties. 

•  Particles  can  be  given  elastic  or  viscoelastic  material  properties;  viscoelastic  material 
properties  are  given  as  a  three-parameter  Maxwell  model. 

•  Cavity  walls  are  assumed  to  be  flat  and  rigid. 

•  A  single  coefficient  of  friction  is  used  for  particle-particle  and  particle-cavity  contact. 

•  The  kinetic  coefficient  of  friction  is  given;  i.e.,  it  is  assumed  that  particle  motion  occurs 
or  that  the  static  coefficient  of  friction  is  zero. 

•  Normal  force-displacement  relations  are  based  on  a  modified  form  of  Lee  &  Radok’s 
relation  (where  the  elastic  modulus  is  replaced  with  the  relaxation  modulus)  with  the 
force  set  to  zero  when  a  negative  force  is  predicted. 

•  For  elastic  properties,  the  normal  force-displacement  relation  reverts  to  an  incremental 
form  of  Hertz’s  relation. 

•  For  viscoelastic  properties,  an  incremental  form  of  the  force-displacement  relation  is 
used  which  incorporates  relaxation  behavior  of  the  viscoelastic  material. 

•  Shear  force-displacement  relations  are  based  on  Amonton’s  law  of  sliding  friction,  or 
Coulomb  friction. 

•  The  shear  force  magnitude  is  based  solely  on  the  magnitude  of  the  normal  force  and  the 
coefficient  of  friction. 

•  The  direction  of  the  shear  force  opposes  the  relative  tangential  velocity  between  the 
contacting  surfaces;  relative  tangential  velocities  can  result  from  oblique  impacts  or  due 
to  rotation  of  the  particles. 

Particle  dampers  are  highly  nonlinear  dampers  whose  energy  dissipation,  or  damping,  is 
derived  from  a  combination  of  loss  mechanisms.  The  relative  effectiveness  of  these 
mechanisms  changes  based  on  various  system  parameters.  The  particle  damper  model  captures 
some  of  the  complex  interactions  involved  in  the  particle  damper.  However,  the  model  tends  to 
overpredict  the  damping.  Results  from  experimental  testing  and  analytical  simulations  of  a 
cantilevered  beam  correlate  reasonably  well.  Correlation  between  experimental  testing  and 
analytical  simulations  of  an  aluminum  chassis  structure  is  somewhat  less.  It  is  believed  that  the 
discrepancy  between  the  experimental  testing  and  analytical  simulations  may  result  largely  due 
to  the  momentum  transfer  that  occurs  when  the  particles  impact  the  cavity  walls.  Any 
compliance  of  the  cavity  walls,  which  are  currently  assumed  to  be  rigid,  may  contribute  to  the 
discrepancy.  It  is  recommended  that  future  efforts  focus  on  ensuring  that  appropriate  forces  are 
applied  during  impacts  between  the  particles  and  the  cavity  walls. 
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APPENDIX  A 

Sample  Input  Deck 


TITLE 

(50)  0.0625"  diameter  steel  particles  with  viscoelastic  properties 
cavity  interior  dimensions:  0.25  x  0.25  x  0.25" 
coefficient  of  friction  =  0.30 
END 

PARAMETERS 

$  TIME  TMAX  TREST  DTMIN  DTMAX  INCR  INCMAX  IREST  NINTPL  ISTAT  ITRAC 
0.  15.000  0.100  5.e-7  5.e-7  0  l.e9  050  2000 
END 

MATID 

$ 

$  Damped  Spring- to -Ground 

$ 

2,  1,  l.Oe-12,  1.0,  1.14882el,  1.00el2,  1.14882el,  1.00el2,  3.08293e-4,  & 
0,  0,  1,  0,  0,  1 
END 

PDMP 

$ 

$  Spherical  Steel  Particle 

$ 

25,  0.03125,  7.32971e-4,  0.33,  25.0e6,  10.0e6,  2.0e-6,  0.30,  0.0,  0.0,  0.0 
END 

BODY 

0.0,  -386.100,  0.0 
END 


TRACE 


25, 

26, 

27, 

28, 

29, 

30, 

31, 

32, 

33, 

34 

Sc 

35,  36, 

37, 

38, 

39, 

40, 

41, 

42, 

43, 

44 

Sc 

45, 

46, 

47, 

48, 

49, 

50, 

51, 

52, 

53, 

54 

Sc 

55, 

56, 

57, 

58, 

59, 

60, 

61, 

62, 

63, 

64 

Sc 

65, 

66, 

67, 

68, 

69, 

70, 

71, 

72, 

73, 

74 

Sc 

75 

END 

NODE  COORDINATES 

$ 

$  Nodes  1-24  define  cavity 
$  Node  25  is  master  cavity  node 
$  Nodes  26-75  define  particle  centers 
$ 

1  O.OOOOOe+000  1.87500e-001  -1 . 25000e-001 

2  O.OOOOOe+000  1.87500e-001  1.25000e-001 

3  2.50000e-001  1.87500e-001  1.25000e-001 

4  2.50000e-001  1.87500e-001  -1 . 25000e-001 

5  O.OOOOOe+000  6.25000e-002  -1 . 25000e-001 

6  O.OOOOOe+000  6.25000e-002  1.25000e-001 

7  2.50000e-001  6.25000e-002  1.25000e-001 

8  2.50000e-001  6.25000e-002  -1 . 25000e-001 

9  1.87500e-001  O.OOOOOe+000  -1 . 25000e-001 

10  1.87500e-001  2.50000e-001  -1 . 25000e-001 

11  1.87500e-001  2.50000e-001  1.25000e-001 

12  1.87500e-001  O.OOOOOe+000  1.25000e-001 

13  6.25000e-002  O.OOOOOe+000  -1 . 25000e-001 

14  6.25000e-002  2.50000e-001  -1 . 25000e-001 

15  6.25000e-002  2.50000e-001  1.25000e-001 

16  6.25000e-002  O.OOOOOe+000  1.25000e-001 

17  O.OOOOOe+000  O.OOOOOe+000  -9 . 37500e-002 

18  2.50000e-001  O.OOOOOe+000  -9 . 37500e-002 

19  2.50000e-001  2.50000e-001  -9 . 37500e-002 

20  O.OOOOOe+000  2.50000e-001  -9 . 37500e-002 

21  O.OOOOOe+000  O.OOOOOe+000  9.37500e-002 
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Best  Available  Copy 


22  2.50000e-001 

23  2.50000e-001 

24  O.OOOOOe+000 

25  O.OOOOOe+000 

26  3.12500e-002 

27  9.37500e-002 

28  1.56250e-001 

29  2.18750e-001 

30  3.12500e-002 

31  9.37500e-002 

32  1.56250e-001 

33  2.18750e-001 

34  3.12500e-002 

35  9.37500e-002 

36  1.56250e-001 

37  2.18750e-001 

38  3.12500e-002 

39  9.37500e-002 

40  1.56250e-001 

41  2.18750e-001 

42  3.12500e-002 

43  9.37500e-002 

44  1.56250e-001 

45  2.18750e-001 

46  3.12500e-002 

47  9.37500e-002 

48  1.56250e-001 

49  2.18750e-001 

50  3.12500e-002 

51  9.37500e-002 

52  1.56250e-001 

53  2.18750e-001 

54  3.12500e-002 

55  9.37500e-002 

56  1.56250e-001 

57  2.187506-001 

58  3.125006-002 

59  9.375006-002 

60  1.562506-001 

61  2.187506-001 

62  3.125006-002 

63  9.375006-002 

64  1.562506-001 

65  2.187506-001 

66  3.125006-002 

67  9.375006-002 

68  1.562506-001 

69  2.187506-001 

70  3.125006-002 

71  9.375006-002 

72  1.562506-001 

73  2.187506-001 

74  3.125006-002 

75  9.375006-002 
END 


0.000006+000 

2.500006-001 

2.500006-001 

0.000006+000 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

9.375006-002 

9.375006-002 

9.375006-002 

9.375006-002 

1.562506-001 

1.562506-001 

1.562506-001 

1.562506-001 

2.187506-001 

2.187506-001 

2.187506-001 

2.187506-001 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

9.375006- 002 

9.375006- 002 

9.375006- 002 

9.375006-002 

1.562506-001 

1.562506-001 

1.562506-001 

1.562506-001 

2.187506-001 

2.187506-001 

2.187506-001 

2.187506-001 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

9.375006-002 

9.375006- 002 

9.375006- 002 

9.375006- 002 

1.562506-001 

1.562506-001 

1.562506-001 

1.562506-001 

2.187506-001 

2.187506-001 

2.187506-001 

2.187506-001 

3.125006-002 

3.125006-002 


9.375006-002 

9.375006-002 

9.375006-002 

O.OOOOOe+000 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9,375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-9.375006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

-3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006- 002 

3.125006- 002 

3.125006- 002 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

3.125006-002 

9.375006-002 

9.375006-002 


ELEMENTS 
1  2  25 
END 


MASS 

25,  3.002E-3,  1.0E20,  1.0E20,  1.0E20 
END 


REEL 

25,123456,  24,  & 

I,  2,  3,  4,  5,  6,  7,  8,  9,  10,  & 

II,  12,  13,  14,  15,  16,  17,  18,  19,  20,  & 
21,  22,  23,  24 

END 
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BOUNDARY  CONDITIONS 
1,  25.  1.  2.  3,  4,  5.  6 
END 

CONTACT 
112  4 
12  3  4 

15  8  6 

16  8  7 

1  9  10  12 
1  10  11  12 
1  13  16  14 
1  14  16  15 
1  17  18  20 
1  18  19  20 
1  21  24  22 
1  22  24  23 
END 

1  26  27  28  29  30  31  32  & 


33 

34 

35 

36 

37 

38 

39 

40 

Sc 

41 

42 

43 

44 

45 

46 

47 

48 

Sc 

49 

50 

51 

52 

53 

54 

55 

56 

Sc 

57 

58 

59 

60 

61 

62 

63 

64 

Sc 

65 

66 

67 

68 

69 

70 

71 

72 

Sc 

73 

74 

75 

END 


Appendix  B 

PDMP  Card 

The  following  pages  give  a  detailed  description  of  the  PDMP  input  block  in  a  format 
similar  to  that  found  in  the  X3D  User’s  Manual  for  the  standard  X3D  input  blocks. 
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PDMP 

INPUT  DATA  BLOCK 


Defines  particle  geometry  and  material  properties  for  particle  damper  simulations. 


STATUS:  Optional  (required  for  particle  damper  simulation) 

FORMAT: 

PDMP 

NODLIM,  RZERO,  PRHO,  PPOI,  EO,  E1 ,  T1 ,  PCOF,  VXO,  VYO,  VZO 
END 


EXAMPLE: 

PDMP 

25,  0.03125,  7.32971  e-4,  0.33,  25.0e6, 10.0e6,  2.0e-6,  0.30,  & 
0.0,  0.0,  0.0 
END 


VARIABLES: 

NODLIM 

RZERO 

PRHO 

PPOI 

EO 

El 

T1 

PCOF 

ViO 


Number  of  last  node  point  in  structure  mesh 

Radius  of  particles 

Density  of  particles 

Poisson’s  ratio  of  particles 

Particle  Eo 

Particle  Ei 

Particle  ti 

Coefficient  of  friction  for  particles 

Initial  velocity  components  for  all  particle  nodes 


NOTES: 

1.  The  particle  nodes  consist  of  nodes  NODLIM+1  through  the  highest 
numbered  node  in  the  model.  All  nodes  in  this  range  are  controlled  by  a 
separate  particle  internal  force  model,  and  should  not  be  connected  to 
other  finite  elements. 
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PDMP 

INPUT  DATA  BLOCK  (CONTINUED) 


2.  RZERO  and  PRHO  define  the  radius  and  density  of  each  particle.  The 
mass  of  a  single  particle  can  be  found  as: 

Mparticle  =  PRHO  X  (4/3)  X  71  X  (RZERO)=^ 

Note  that  if  the  highest  numbered  node  in  the  model  is  Nmax,  then  the  total 
mass  of  the  particles  can  be  found  as: 


Mtotal  =  [Nmax  ~  NODLIM]  X  Mparticle 

3.  PPOl,  EO,  El ,  and  T1  define  the  material  properties  of  the  particles.  PPOl 
is  set  to  the  Poisson’s  ratio  for  the  material.  For  elastic  properties,  EO  is  set 
to  the  elastic  modulus,  E1  is  set  to  zero,  and  T 1  is  set  to  any  nonzero 
number  (i.e.,  Tt^^O).  For  viscoelastic  properties,  EO,  E1,  and  T1  define  the 
parameters  in  a  three-parameter  Maxwell  model  representation. 

4.  PCOF  defines  the  coefficient  of  friction.  A  single  coefficient  of  friction  is 
used  for  particle-particle  and  particle-cavity  contact.  The  kinetic  coefficient 
of  friction  is  given;  i.e.,  it  is  assumed  that  particle  motion  occurs  or  that  the 
static  coefficient  of  friction  is  zero. 

5.  VXO,  VYO,  and  VZO  define  the  initial  components  of  velocity  for  all 
particle  nodes.  Nonuniform  initial  velocity  values  may  be  input  using  the 
INITial  input  block.  All  initial  velocities  are  subject  to  the  boundary 
conditions  prescribed  in  the  BOUN  input  block. 
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Appendix  C 

Sample  USRFRC  Routine 


Q  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  ** 

Q  **  ★*  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  **  USRFRC 

C 

SUBROUTINE  USRFRC  (  NUMNOD,  ncord/  XYZINI,  XYZ,  VBC,  VEL,  FORCE  ) 

C 

C  Hook  for  user- supplied  subroutine  to  apply  nodal  forces  at  each  time 
C  step.  Called  from  CNTDIF  at  start  of  solution  step.  The  force  array 
C  contains  (external  force  -  internal  force)  for  each  DOF. 

C 

C  *  NUMNOD  =  Number  of  nodes  in  model  (DO  NOT  MODIFY) 

C  *  ncord  =  Leading  dimension  of  coordinate  arrays  (3  or  6) 

C  *  XYZINI  =  Initial  nodal  coordinates  (DO  NOT  MODIFY) 

C  *  XYZ  =  Coordinates  at  last  time  step  (DO  NOT  MODIFY) 

C  *  VBC  =  Velocity  BC  codes  (DO  NOT  MODIFY) 

C  *  VEL  =  Current  velocity  values  (DO  NOT  MODIFY) 

C  *  FORCE  =  Current  force  sums  for  each  nodal  DOF 
C 

C  To  apply  a  prescribed  nodal  force,  ADD  the  corresponding  force  to  the 
C  nodal  force  array. 

C 

C  WARNING:  Data  are  not  stored  in  nodal  arrays  by  * external* 

C  node  number,  but  by  sequential  "internal"  numbers  which  are 
C  invisible  to  the  user. 

C 

C  To  get  the  internal  node  number  for  a  node  with  user-assigned 
C  node  number  "NID",  use  CALL  GETIDN  (  NID,  IDSEQ  ) .  The  n\mnber 
C  IDSEQ  can  be  used  to  address  nodal  data  in  arrays. 

C 

C  Example:  to  assign  an  X-force  of  50  units  to  node  545: 

C 

C  NID  =  545 

C  CALL  GETIDN  (  NID,  IDSEQ  ) 

C  FORCE (1, IDSEQ)  =  FORCE ( 1 , IDSEQ)  +  50. 

C 

C  Forces  are  positive  in  the  positive  coordinate  direction. 

C 

include  'x3d.ins' 
c 

COMMON  /  SOLPAR  /  TIME  ,  TMAX  ,  TREST  ,  DTMIN  ,  DTMAX  , 

+  INCR  ,  INCMAX,  IREST  ,  ITRACE,  IMPOST 
c 

DIMENSION  XYZINI (ncord, *) ,  XYZ (ncord, *) ,  VEL(6,*),  FORCE (6,*) 
CHARACTER*!  VBC (6,*) 
c 

c  37.21  Hz  =  233.7973253  rad/sec 
c 

FRCNEW  =  4.84040e-3  *  SIN  (233.7973253  *  TIME) 
c 

CALL  GETIDN  (  25,  NODE  ) 

FORCE (1, NODE)  =  FORCE ( 1 , NODE )  +  FRCNEW 
c 

RETURN 

END 
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Appendix  D: 

MATLAB  Cavity  Fill  Program 


%  create_cavity.in  --  Create  SDOF  cavity-particle  geometry 
%  01-Dec -2 001  CSA  Engineering 

%  Note:  stacking  algorithm  is  simplistic  and  does  not  calculate  packing 
%  correctly.  It  works  for  most  cases  of  spherical  particles,  though, 
clear  all; 

%  Enter  cavity  dimensions 
x=0.25;  y=0.25;  z=^0.25; 

%  Enter  particle  diameter 
diam  =  0.0625; 

%  Enter  number  of  particles 
n  =  50; 

% 

%  Do  not  change  below  here 
% 

px  =  zeros {n,l); 
py  =  zeros (n, 1) ; 
pz  =  zeros (n,l); 


cav 

=  [... 

1 

0  3*y/4 

-z/2  ; 

2 

0  3*y/4 

z/2  ;. 

3 

X  3*y/4 

z/2  ;. 

4 

K  3*y/4 

-z/2  ; 

5 

0  y/4  - 

z/2  ;., 

6 

0  y/4  z 

/2  ;... 

7 

X  y/4  z 

/2  ;... 

8 

X  y/4  - 

z/2  ;.. 

9 

3*x/4  0 

-z/2  ; 

10 

3*x/4 

y  -z/2 

11 

3*x/4 

y  z/2  ; 

12 

3*x/4 

0  z/2  ; 

13 

x/4  0 

-z/2  ;  . 

14 

x/4  y 

-z/2  ;  . 

15 

x/4  y 

z/2  ;.. 

16 

x/4  0 

z/2  ;.. 

17 

0  0-3 

*z/8;  .  . 

18 

X  0  -3 

*z/8;  .  . 

19 

X  y  -3 

*z/8;  .  . 

20 

0  y  -3 

*z/8;  .  . 

21 

0  0  3* 

z  /  8 ;  .  .  . 

22 

X  0  3* 

z  /  8 ;  .  .  . 

23 

X  y  3* 

z/8;  .  .  . 

24 

0  y  3* 

z  /  8 ;  .  .  . 

25 

0  0  0] 

f 

r  =  diam/2; 
px(l:n)  =  r; 
py(l:n)  =  r; 
pz (1 :n)  =  -z/2  +  r ; 
for  indx  =  2:n 

px(indx)  =  px(indx-l)  +  diam; 
if  (px(indx)  >  x) 
px(indx:n)  =  r; 

py{indx:n)  =  py(indx-l)  +  diam; 
if  (py(indx)  >  y) 
px{indx:n)  =  r; 
py(indx:n)  =  r; 

pz(indx:n)  =  pz(indx-l)  +  diam; 
if  (pz(indx)  >  z) 

disp ([ 'Error :  particle  '  num2str (indx)  '  outside  cavity.']); 
disp( [ 'Reduce  n  to  something  less  than  this  and  rerun.']); 
return; 
end 
end 
end 
end 

index  =  26:n+25; 
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index  =  index ' ; 
p  =  [  index,  px,  py,  pz  ] ; 

xc=cav( : ,2) ; 
yc^cav { : , 3 ) ; 
zc=cav{ : , 4) ; 
xp=p ( : , 2 ) ; 
yp=p ( : / 3 )  ; 
zp=p ( : , 4 ) ; 

plot3  (xc,yc,zc,  'k+'  ,xp,yp,zp,  ; 

grid  on; 

xlabel ( 'x' ) ; 

ylabel ( 'y' ) ; 

zlabel ( ' z ' ) ; 

fprintf ( ' \nCopy/Paste  the  following  table  into  your  x3d  NODES. \nM; 
for  indx  =  1:25 

fprintf(  '  %d  %10.5e  %10.5e  %10 . 5e\nS cav(indx, : ) ) ; 
end 

for  indx  =  l:n 

fprintf(  '  %d  %10.5e  %10.5e  %10 . 5e\n^p  (indx,  : )  )  ; 
end% 
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Appendix  E 

MATLAB  Post  Processing  Program 


% 

%  pmov.m 
% 

%  MATLAB  routine  to  make  particle  damper  movies 
%  using  results  from  X3D  trace  files. 

% 

%  11  Dec  98  -  SEO  @  UDRI 

% 

cir=0:pi/40:2*pi;  %  used  to  draw  spheres 

nump  =64;  %  number  of  particles 
prad  =  0.03125;  %  particle  radius 
%clen  =  1.00;  %  cavity  length 
%chgt  =  1.00;  %  cavity  height 
clen  =  0.263;  %  cavity  length 
chgt  =  0.375;  %  cavity  height 
cwdt  =  0.375;  %  cavity  width 

% 

%  open  trace  file  and  read  data 
% 


fid  =  fopen( 'trace.txt' , 'r' ) ; 
temp  =  dlmread { ' trace . txt 

% 

%  calculate  number  of  rows  and  times  where  data 
%  is  available 
% 

rows  =  si2e(temp, 1) ;  %  number  of  row  in  trace  file 
tims  =  rows/ (nump+1) ;  %  number  of  times  at 
%  which  data  is  avialable 

maxv  =  0.0;  %  set  max  trans  velocity  to  0.0 
maxr  =  0.0;  %  set  max  rot  velocity  to  0.0 

% 

%  process  data  as  read  from  trace  file 
% 


tstep  =  20; 

for  j=l: tstep: tims 
simt{j)  =  temp ( ( j-1)  *  (nuitrp+1) +1, 1)  ; 
for  k=l:nump+l 

% 

%  get  particles'  position  (x,y,z) 

% 

ppos ( j ,k, 1) =temp( ( j-1) * (nump+1) +k/ 2} ; 
ppos ( j ,k, 2) =temp( (j-1) * (nump+1) +k/ 3) ; 
ppos ( j ,k, 3) =temp( (j-1) * (nump+1) +k, 4) ; 

% 

%  get  particles'  translational  velocity 
% 

pvel  ( j  ;kwl)  =temp{  (j-1)  *  (nump+1) +k,  5) ; 
pvel(j/k/2)=temp(  (j-l)  *  (nuirp+l) +k,  6)  ; 

% 

%  calculate  maximum  translational  velocity 
% 

mvel=sqrt (pvel ( j , k/ 1) ^2+pvel { j ,k/ 2) ^2) ; 
if  mvel  >  maxv 
maxv=mvel ; 
end 
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% 

%  get  particles'  rotational  velocity 
% 

pvel  { j  /  k,  3)=teTnp(  ( j-1)  *  (nump+l)  +k,  7)  ; 

% 

%  calculate  maximum  rotational  velocity 
% 

if  abs(pvel{j,k,3) )  >  maxr 
maxr=abs (pvel ( j / k^  3 ) ) ; 
end 
end 
end 


% 

%  set  axes  for  plotting  simulation  data 
%  and  set  aside  space  for  movie 
% 

axis  equal; 

%axis([-1.00  1.00  -1.00  1.00]); 
axis([-0.20  0.50  -0.10  0.60]); 

%axis([-1.00  2.00  -1.00  4.00]); 

% 

%  moviein  no  longer  required  in  Matlab  5 . 3 
% 

%  M  =  moviein  ( tims ) ; 

set (gca, 'NextPlot' , 'replacechildren' ) 

% 

%  loop  over  the  times  at  which  data  available 
% 

for  j=l : tstep: tims 
% 

%  clear  existing  figure 
% 

cla 

% 

%  draw  cavity's  exterior  walls 
% 

%  x=[ppos ( j , 1, 1) -0.100;ppos( j ,1,1) +clen+0.100] ; 

%  x= [x;ppos (j,l,l)+clen+0. 100 ;ppos(j, 1,1) -0.100] ; 

%  x= [x;ppos ( j , 1, 1) -0. 100] ; 

%  y= [ppos { j , 1, 2 ) -0 . 100;ppos ( j , 1 , 2 ) -0 . 100] ; 

%  y= [y;PPOS { j , l , 2 ) +chgt+0 . 100 ;ppos { j , 1 , 2 ) +chgt+0 . 100 ] ; 

%  y= [y;ppos { j , 1 , 2 ) -0 . 100] ; 
x= [ppos ( j , 1 , 1) -0 . 062 ;ppos ( j , 1 , 1) +clen+0 . 062] ; 
x= [x; ppos (j ,1,1) +clen+0. 062 ;ppos (j,l,l) -0.062] ; 
x= [x;ppos ( j , 1, 1) -0 . 062] ; 
y= [ppos { j , 1 , 2) -0 . 064 ; ppos {j , 1 , 2) -0 . 064] ; 
y= [y;PPOS ( j , l , 2 ) +chgt+0 . 064 ;ppos ( j , 1 , 2) +chgt+0 . 064] ; 
y= [y;ppos ( j , 1, 2) -0 . 064] ; 

% 

%  x= [ppos { j , 1, 1) ;ppos ( j , 1, l)+clen+0.124] ; 

%  x=[x;ppos{j,l,l)+clen+0.124;ppos(j,l,l) ;ppos(j,l,l) ] 
%  y= [ppos (j, 1,2) ;ppos(j,l,2) ;ppos(j,l,2)+chgt+0.128] ; 

%  y= [y;ppos( j ,1,2) +chgt+0.128;ppos{ j ,1, 2) ] ; 

%  plot(x,y) 
patch(x,y, [0  0  1] ) 
hold  on 

% 

%  draw  cavity's  interior  walls 
% 


x= [ppos ( j , 1 , 1 ) ;ppos ( j , 1 , 1 ) +clen ;ppos { j , 1 , 1 ) +clen] ; 
x= [x;ppos ( j , 1, 1) ;ppos {j,l,l)]; 
y= [ppos (j, 1,2) ;ppos(j,l,2) ;ppos ( j , 1 , 2 ) +chgt] ; 
y=[y;ppos(j,l,2)+chgt;ppos{j,l,2) ] ; 

%  x= [ppos { j, l,l)+0. 062 ;ppos(j,l,l)+clen+0. 062] ; 

%  x= [x;ppos { j, 1, 1 ) +clen+0. 062 ;ppos(j, 1,1) +0.062] ; 
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%  x=[x;ppOs(j ,l,l)+0.062] ; 

%  y=[ppos(j,l/2)+0.064;ppos{j,l,2)+0.064] ; 

%  y= [y ;PPOS {3,1,2) +chgt+0 . 064 ;ppos ( j , 1 , 2 ) +chgt+0 . 064 ] ; 

%  y= [y;ppos ( j , 1, 2) +0 . 064] ; 

%  plot(x,y) 
patch(X/y, 'w' ) 
axis  equal; 

%  axis([-1.00  1.00  -1.00  1.00]); 

axis(E-0.20  0.50  -0.10  0.60]); 

%  axis{[-1.00  2.00  -1.00  4.00]); 

% 

%  sort  particles  in  z-direction 
% 

for  k=2:nuiTp 

for  l=k+l  :nuitp+l 

if  ppos(j,lr3)  <  ppos(j/k,3) 

tempi  =  ppos { j / 1, 1) ; 

temp2  =  ppos ( j ,1,2); 

temp3  =  ppos ( j , 1,3) ; 

ppos{j,l,l)  =  ppos(j,k,l) ; 

ppos(j,l,2)  =ppos(j,k,2); 

ppos (j, 1,3)  =  ppos(j,k,3); 

ppos(j,k,l)  =  tempi; 

ppos(j,k,2)  =  temp2; 

ppos ( j , k, 3 )  =  temp3 ; 

end 

end 

end 

% 

%  draw  particle  boundaries 
% 

for  k=2:nump+l 

a=[ppos ( j ,k, l)+prad*sin(cir) ] ; 
b=Eppos(j,k,2)+prad*cos(cir) ] ; 

% 

%  for  filled  particles 

%  (currently  set  up  for  four  yellow  particles) 

% 

val  -  ppos(j,k,3)  +  (  cwdt  /  2  ); 
val  =  val  /  (  cwdt  ) ; 

% 

%  if  k=:=2 

patch(a,b, [1  val  val]) 

%  end 
%  if  k==3 
%  patch{a,b, 'g' ) 

%  end 
%  if  k==4 
%  patch(a,b,  'y' ) 

%  end 
%  if  k==5 
%  patch{a,b, 'y' ) 

%  end 
% 

%  for  unfilled  particles 
% 

%  plot (a, b) 

% 

%  draw  point  at  particle  center 
% 

%  plot (ppos (j ,k, 1) ,ppos(j,k,2) , ' 

%  .  _  . 

%  draw  line  from  particle  center  to  outside  to  track  rotation 

% 

%  a=[ppos(j,k,l)  ppos(j,k,l)+prad*cos(ppos(j,k,3) ) ] ; 

%  b=[ppos(j,k,2)  ppos ( j ,k, 2)+prad*sin(ppos (j ,k, 3) ) ] ; 

% 

%  plot (a,b, 'b' ) 

% 

%  draw  line  to  show  translational  velocity  at  current  time 
% 
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%  vmag  =  sqrt  {pvel  ( j  1)  ^2+pvel  ( j /k,  2) '^2)  ; 

% 

%  if  (vmag/maxv  >  0.20) 

%  sfac  =  prad  *  (vinag  /  maxv)  ; 

%  xmag  =  sfac  *  pvel  { j  ,  k,  1) /vmag; 

%  ymag  =  sfac  *  pvel { j , k, 2) /vmag; 

%  scale  =  1.5; 

%  theta  =  at an (ymag/xmag) ; 

%  if  (xmag  >=  0) 

%  scale  =  -scale; 

%  end 

%  a= [ppos ( j ,k/ 1)  ppos ( j /k/l)+xmag] ; 

%  a=[a  ppos{j,k,l)+xmag+0.01*scale*cos(theta+pi/6) ] ; 

%  a=[a  ppos ( j / k/ 1) +xmag3 ; 

%  a=[a  ppos{j,k,l)+xmag+0.01*scale*cos(theta-pi/6} ] ; 

%  b=[ppos(j,k,2)  ppos(j,k,2)+ymag] ; 

%  b=[b  ppos(j,k,2)+ymag+0.01*scale*sin(theta+pi/6) ] ; 

%  b=[b  ppos ( j ,k,2)+ymag] ; 

%  b=[b  ppos(j,k,2)+ymag+0.01*scale*sin(theta-pi/6) ] ; 

% 

%  plot (a/b/ 'r' ) 

% 

%  end 

draw  arc  to  show  rotational  velocity  at  current  time 

%  angle  =  pi* (pvel { j /k, 3) /maxr) ; 

%  if  (abs (angle)  >  pi/9) 

%  arc  =  0 : angle/20: angle; 

%  xstr  =  ppos(j,k,l); 

%  ystr  =  ppos(j,k,2); 

%  scale  =  1.5; 

%  if  (pvel(j,k,3)  >=  0) 

%  theta=angle-pi/2+pi/18; 

%  a=:  [xstr+0.5*prad*cos  (arc)  ]  ; 

%  a=[a  xstr+0. 5*prad*cos (angle) +0.01*scale*cos (theta) 1 ; 

%  b= [ystr+0 . 5*prad*sin(arc) ] ; 

%  b=[b  ystr+0. 5*prad*sin(angle)+0.01*scale*sin(theta) ] ; 

%  end 

%  if  (pvel(j,k,3)  <  0) 

%  theta=angle+pi/2-pi/18; 

%  a= [xstr+0. 5*prad*cos (arc) ] ; 

%  a=[a  xstr+0. 5*prad*cos (angle) +0.01*scale*cos (theta) ] ; 

%  b= [ystr-0 . 5*prad*sin{arc) ] ; 

%  b=[b  ystr-0.5*prad*sin(angle)+0,01*scale*sin(theta) ] ; 

%  end 

plot (a/b/ 'r') 

%  end 
end 

label  plot  with  simulated  time 

labt  =  [  'Simulated  Time  (sec):  num2str (simt  ( j ) / ' %0 , 5e' )  ] 

%  text (  -0.800,  -0.800/  labt) 
text(  -0.15,  -0.05,  labt) 

%  text(  -0.50/  -0.50/  labt) 

labt  =  'Particle  Damper  Simulation'; 

%  text{  -0.800,  0.800/  labt) 
text (  -0.15/  0.55/  labt) 

%  text(  -0.50/  -0.75/  labt) 
hold  off 

m  =  int8 ( ( ( j -1) /tstep) +1) ; 

M(:,m)  =  getframe; 

P  =  getframe; 

%  tname  =  ['frame'  int2str (m) ] ' ; 

%  imwrite (P.cdat a, tname, 'bmp' ) ; 
if  j==l 
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fname  =  ^frOOOO^* 
print ( ' -dmeta ' , fname ) ; 

%  print (' -dmfile' , fname) ; 

%  imwri te ( P . cdata , fname ,  ' emf ' ) ; 
end 
if  j>l 

if  j/tstep  <  1000 

fname  =  ['frO'  int2str( j/tstep) ] ; 

end 

if  j/tstep  <  100 

fname  =  ['frOO'  int2str( j/tstep) ] ; 
end 

if  j/tstep  <  10 

fname  =  ['frOOO'  int2str { j /tstep) ] ; 
end 

print ( *  -dmeta' , fname) ; 

% 

%  print  (' -dmfile' ,  fname)  ; 

%  imwri te (P. cdata, fname, 'emf M ; 
end 

%  if  rem(j“l,5)  ==  0.0 
%  fname  =  ['frame'  int2str ( j-1) ] ; 

%  print (' -dmfile ', fname) ; 

%  end 

end 
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List  of  Acronyms 


ACRONYM 

DESCRIPTION 

AFRL 

Air  Force  Research  Laboratory 

FE 

finite  element 

FFT 

fast  fourier  transform 

FRF 

frequency  response  function 

MDOF 

multiple-degree-of-freedom 

MIP 

multiple  individual  particle 

MPDD 

multiparticle  impact  damping 

MSD 

mass-spring-damper 

PDMP 

X3D  particle  damper  card 

SBIR 

Small  Business  Innovation  Research 

SDOF 

single-degree-of-freedom 

STIR 

Small  Business  Technology  Transfer 

UDRI 

University  of  Dayton  Research  Institute 

USAF 

United  States  Air  Force 
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